CellStepper.C 1.97 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2011 The Regents of the University of California
//
// This file is part of Qbox
//
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// CellStepper.C
//
////////////////////////////////////////////////////////////////////////////////

#include "CellStepper.h"
using namespace std;

////////////////////////////////////////////////////////////////////////////////
23
void CellStepper::enforce_constraints(double *u)
24
{
25
  // modify the strain u so as to satisfy the constraints expressed in cell_lock
26 27 28 29 30 31 32 33
  string cell_lock = s_.ctrl.cell_lock;

  if ( cell_lock != "OFF" )
  {
    // constraints on the cell derivatives
    if ( cell_lock.find("A") != string::npos )
    {
      // vector A is locked
34
      u[0] = 0.0;
35 36 37 38
    }
    if ( cell_lock.find("B") != string::npos )
    {
      // vector B is locked
39
      u[4] = 0.0;
40 41 42 43
    }
    if ( cell_lock.find("C") != string::npos )
    {
      // vector C is locked
44
      u[8] = 0.0;
45 46 47 48 49 50
    }

    // Check if cell shape should be preserved (if "S" is present in cell_lock)
    // The only changes allowed are renormalizations of a,b,c
    if ( cell_lock.find("S") != string::npos )
    {
51 52
      // set all off-diagonal elements of u to zero
      u[1] = u[2] = u[3] = u[5] = u[6] = u[7] = 0.0;
53 54 55 56
    }

    if ( cell_lock == "R" )
    {
57 58 59 60 61
      // set all off-diagonal elements of u to zero
      u[1] = u[2] = u[3] = u[5] = u[6] = u[7] = 0.0;
      // set all diagonal elements to 1/3 trace u
      const double tr3rd = ( u[0] + u[4] + u[8] ) / 3.0;
      u[0] = u[4] = u[8] = tr3rd;
62 63 64
    }
  }

65
  // u is a strain tensor respecting constraints
66
}