v0.14.0
InterfaceGapArcLengthControl.hpp
Go to the documentation of this file.
1 /** \file InterfaceGapArcLengthControl.hpp
2  \brief Implementation of arc-lebgth control for cohesive elements
3 
4  Arc-length in that version controls gap opening
5 
6 */
7 
8 
9 
10 namespace CohesiveElement {
11 
12 struct ArcLengthIntElemFEMethod: public FEMethod {
13  moab::Interface& mOab;
14 
15  boost::shared_ptr<ArcLengthCtx> arcPtr;
16 
17  Vec GhostLambdaInt;
18  Range Faces3,Faces4;
19  Range Edges3,Edges4;
20  Range Nodes3,Nodes4;
21 
22  Tag thDamagedPrism;
23 
24  ArcLengthIntElemFEMethod(
25  moab::Interface& moab,
26  boost::shared_ptr<ArcLengthCtx> arcptr):
27  FEMethod(),mOab(moab),arcPtr(arcptr) {
28  PetscInt ghosts[1] = { 0 };
29  ParallelComm* pcomm = ParallelComm::get_pcomm(&mOab,MYPCOMM_INDEX);
30  if(pcomm->rank() == 0) {
31  VecCreateGhost(PETSC_COMM_WORLD,1,1,0,ghosts,&GhostLambdaInt);
32  } else {
33  VecCreateGhost(PETSC_COMM_WORLD,0,1,1,ghosts,&GhostLambdaInt);
34  }
35  Range prisms;
36  rval = mOab.get_entities_by_type(0,MBPRISM,prisms,false); MOAB_THROW(rval);
37  for(Range::iterator pit = prisms.begin();pit!=prisms.end();pit++) {
38  EntityHandle f3,f4;
39  rval = mOab.side_element(*pit,2,3,f3); MOAB_THROW(rval);
40  rval = mOab.side_element(*pit,2,4,f4); MOAB_THROW(rval);
41  Faces3.insert(f3);
42  Faces4.insert(f4);
43  }
44 
45  rval = mOab.get_adjacencies(Faces3,1,false,Edges3); MOAB_THROW(rval);
46  rval = mOab.get_adjacencies(Faces4,1,false,Edges4); MOAB_THROW(rval);
47  rval = mOab.get_connectivity(Faces3,Nodes3,true); MOAB_THROW(rval);
48  rval = mOab.get_connectivity(Faces4,Nodes4,true); MOAB_THROW(rval);
49  //Faces3.insert(Edges3.begin(),Edges3.end());
50  Faces3.insert(Nodes3.begin(),Nodes3.end());
51  //Faces4.insert(Edges4.begin(),Edges4.end());
52  Faces4.insert(Nodes4.begin(),Nodes4.end());
53 
54  double def_damaged = 0;
55  rval = mOab.tag_get_handle(
56  "DAMAGED_PRISM",1,MB_TYPE_INTEGER,thDamagedPrism,MB_TAG_CREAT|MB_TAG_SPARSE,&def_damaged); MOAB_THROW(rval);
57 
58  }
59  ~ArcLengthIntElemFEMethod() {
60  VecDestroy(&GhostLambdaInt);
61  }
62 
63  /** \brief remove nodes of prims which are fully damaged
64  *
65  */
66  MoFEMErrorCode remove_damaged_prisms_nodes() {
67  MoFEMFunctionBeginHot;
68  Range prisms;
69  rval = mOab.get_entities_by_type(0,MBPRISM,prisms,false); CHKERRG(rval);
70  std::vector<int> is_prism_damaged(prisms.size());
71  rval = mOab.tag_get_data(thDamagedPrism,prisms,&*is_prism_damaged.begin()); CHKERRG(rval);
72  Range::iterator pit = prisms.begin();
73  std::vector<int>::iterator vit = is_prism_damaged.begin();
74  for(;pit!=prisms.end();pit++,vit++) {
75  if(*vit>0) {
76  Range nodes;
77  rval = mOab.get_connectivity(&*pit,1,nodes,true); CHKERRG(rval);
78  for(Range::iterator nit = nodes.begin();nit!=nodes.end();nit++) {
79  Faces3.erase(*nit);
80  Faces4.erase(*nit);
81  }
82  }
83  }
84  MoFEMFunctionReturnHot(0);
85  }
86 
87  double lambda_int;
88  MoFEMErrorCode preProcess() {
89  MoFEMFunctionBeginHot;
90  switch(snes_ctx) {
91  case CTX_SNESSETFUNCTION: {
92  ierr = calculate_dx_and_dlambda(snes_x); CHKERRG(ierr);
93  ierr = calculate_db(); CHKERRG(ierr);
94  ierr = calculate_lambda_int(lambda_int); CHKERRG(ierr);
95  }
96  break;
97  default:
98  break;
99  }
100  MoFEMFunctionReturnHot(0);
101  }
102 
103  MoFEMErrorCode calculate_lambda_int(double &_lambda_int_) {
104  MoFEMFunctionBeginHot;
105  ParallelComm* pcomm = ParallelComm::get_pcomm(&mOab,MYPCOMM_INDEX);
106  NumeredDofEntityByLocalIdx::iterator dit,hi_dit;
107  dit = problemPtr->getNumeredRowDofsPtr()->get<PetscLocalIdx_mi_tag>().lower_bound(0);
108  hi_dit = problemPtr->getNumeredRowDofsPtr()->get<PetscLocalIdx_mi_tag>().upper_bound(problemPtr->getNbLocalDofsRow());
109  double *array;
110  double *array_int_lambda;
111  ierr = VecZeroEntries(GhostLambdaInt); CHKERRG(ierr);
112  ierr = VecGhostUpdateBegin(GhostLambdaInt,INSERT_VALUES,SCATTER_FORWARD); CHKERRG(ierr);
113  ierr = VecGhostUpdateEnd(GhostLambdaInt,INSERT_VALUES,SCATTER_FORWARD); CHKERRG(ierr);
114  ierr = VecGetArray(arcPtr->dx,&array); CHKERRG(ierr);
115  ierr = VecGetArray(GhostLambdaInt,&array_int_lambda); CHKERRG(ierr);
116  array_int_lambda[0] = 0;
117  for(;dit!=hi_dit;dit++) {
118  if(dit->get()->getEntType() != MBVERTEX) continue;
119  if(pcomm->rank() != dit->get()->getPart()) continue;
120  if(Nodes3.find(dit->get()->getEnt())!=Nodes3.end()) {
121  array_int_lambda[0] += array[dit->get()->getPetscLocalDofIdx()];
122  }
123  if(Nodes4.find(dit->get()->getEnt())!=Nodes4.end()) {
124  array_int_lambda[0] -= array[dit->get()->getPetscLocalDofIdx()];
125  }
126  }
127  //PetscSynchronizedPrintf(PETSC_COMM_WORLD,
128  // "array_int_lambda[0] = %6.4ee\n",array_int_lambda[0]);
129  ierr = VecRestoreArray(arcPtr->dx,&array); CHKERRG(ierr);
130  ierr = VecRestoreArray(GhostLambdaInt,&array_int_lambda); CHKERRG(ierr);
131  ierr = VecGhostUpdateBegin(GhostLambdaInt,ADD_VALUES,SCATTER_REVERSE); CHKERRG(ierr);
132  ierr = VecGhostUpdateEnd(GhostLambdaInt,ADD_VALUES,SCATTER_REVERSE); CHKERRG(ierr);
133  ierr = VecGhostUpdateBegin(GhostLambdaInt,INSERT_VALUES,SCATTER_FORWARD); CHKERRG(ierr);
134  ierr = VecGhostUpdateEnd(GhostLambdaInt,INSERT_VALUES,SCATTER_FORWARD); CHKERRG(ierr);
135  ierr = VecGetArray(GhostLambdaInt,&array_int_lambda); CHKERRG(ierr);
136  _lambda_int_ = arcPtr->alpha*array_int_lambda[0] + arcPtr->dLambda*arcPtr->beta*sqrt(arcPtr->F_lambda2);
137  /*PetscSynchronizedPrintf(PETSC_COMM_WORLD,
138  "array_int_lambda[0] = %6.4e arcPtr->F_lambda2 = %6.4e\n",
139  array_int_lambda[0],arcPtr->F_lambda2);
140  PetscSynchronizedFlush(PETSC_COMM_WORLD);*/
141  ierr = VecRestoreArray(GhostLambdaInt,&array_int_lambda); CHKERRG(ierr);
142  MoFEMFunctionReturnHot(0);
143  }
144 
145  virtual MoFEMErrorCode calculate_db() {
146  MoFEMFunctionBeginHot;
147  ierr = VecZeroEntries(arcPtr->db); CHKERRG(ierr);
148  ierr = VecGhostUpdateBegin(arcPtr->db,INSERT_VALUES,SCATTER_FORWARD); CHKERRG(ierr);
149  ierr = VecGhostUpdateEnd(arcPtr->db,INSERT_VALUES,SCATTER_FORWARD); CHKERRG(ierr);
150  NumeredDofEntityByLocalIdx::iterator dit,hi_dit;
151  dit = problemPtr->getNumeredRowDofsPtr()->get<PetscLocalIdx_mi_tag>().lower_bound(0);
152  hi_dit = problemPtr->getNumeredRowDofsPtr()->get<PetscLocalIdx_mi_tag>().upper_bound(
153  problemPtr->getNbLocalDofsRow()+problemPtr->getNbGhostDofsRow()
154  );
155  double *array;
156  ierr = VecGetArray(arcPtr->db,&array); CHKERRG(ierr);
157  for(;dit!=hi_dit;dit++) {
158  if(dit->get()->getEntType() != MBVERTEX) {
159  array[dit->get()->getPetscLocalDofIdx()] = 0;
160  continue;
161  }
162  if(Nodes3.find(dit->get()->getEnt())!=Nodes3.end()) {
163  array[dit->get()->getPetscLocalDofIdx()] = +arcPtr->alpha;
164  } else if(Nodes4.find(dit->get()->getEnt())!=Nodes4.end()) {
165  array[dit->get()->getPetscLocalDofIdx()] = -arcPtr->alpha;
166  } else array[dit->get()->getPetscLocalDofIdx()] = 0;
167  }
168  ierr = VecRestoreArray(arcPtr->db,&array); CHKERRG(ierr);
169  MoFEMFunctionReturnHot(0);
170  }
171 
172  MoFEMErrorCode operator()() {
173  MoFEMFunctionBeginHot;
174 
175  switch(snes_ctx) {
176  case CTX_SNESSETFUNCTION: {
177  //calculate residual for arc length row
178  arcPtr->res_lambda = lambda_int - arcPtr->s;
179  ierr = VecSetValue(snes_f,arcPtr->getPetscGlobalDofIdx(),arcPtr->res_lambda,ADD_VALUES); CHKERRG(ierr);
180  PetscPrintf(PETSC_COMM_SELF,"\tres_lambda = %6.4e lambda_int = %6.4e s = %6.4e\n",
181  arcPtr->res_lambda,lambda_int,arcPtr->s);
182  }
183  break;
184  case CTX_SNESSETJACOBIAN: {
185  //calculate diagonal therm
186  arcPtr->dIag = arcPtr->beta*sqrt(arcPtr->F_lambda2);
187  ierr = MatSetValue(snes_B,arcPtr->getPetscGlobalDofIdx(),arcPtr->getPetscGlobalDofIdx(),1,ADD_VALUES); CHKERRG(ierr);
188  }
189  break;
190  default:
191  break;
192  }
193 
194  MoFEMFunctionReturnHot(0);
195  }
196 
197  MoFEMErrorCode postProcess() {
198  MoFEMFunctionBeginHot;
199  switch(snes_ctx) {
200  case CTX_SNESSETJACOBIAN: {
201  ierr = VecGhostUpdateBegin(arcPtr->ghostDiag,INSERT_VALUES,SCATTER_FORWARD); CHKERRG(ierr);
202  ierr = VecGhostUpdateEnd(arcPtr->ghostDiag,INSERT_VALUES,SCATTER_FORWARD); CHKERRG(ierr);
203  PetscPrintf(PETSC_COMM_WORLD,"\tdiag = %6.4e\n",arcPtr->dIag);
204  ierr = MatAssemblyBegin(snes_B,MAT_FLUSH_ASSEMBLY); CHKERRG(ierr);
205  ierr = MatAssemblyEnd(snes_B,MAT_FLUSH_ASSEMBLY); CHKERRG(ierr);
206  }
207  break;
208  default:
209  break;
210  }
211  MoFEMFunctionReturnHot(0);
212  }
213 
214  MoFEMErrorCode calculate_dx_and_dlambda(Vec &x) {
215  MoFEMFunctionBeginHot;
216  //dx
217  ierr = VecCopy(x,arcPtr->dx); CHKERRG(ierr);
218  ierr = VecAXPY(arcPtr->dx,-1,arcPtr->x0); CHKERRG(ierr);
219  //if LAMBDA dof is on this partition
220  if(arcPtr->getPetscLocalDofIdx()!=-1) {
221  double *array;
222  ierr = VecGetArray(arcPtr->dx,&array); CHKERRG(ierr);
223  arcPtr->dLambda = array[arcPtr->getPetscLocalDofIdx()];
224  array[arcPtr->getPetscLocalDofIdx()] = 0;
225  ierr = VecRestoreArray(arcPtr->dx,&array); CHKERRG(ierr);
226  }
227  //brodcast dlambda
228  ierr = VecGhostUpdateBegin(arcPtr->ghosTdLambda,INSERT_VALUES,SCATTER_FORWARD); CHKERRG(ierr);
229  ierr = VecGhostUpdateEnd(arcPtr->ghosTdLambda,INSERT_VALUES,SCATTER_FORWARD); CHKERRG(ierr);
230  //calculate dx2 (dot product)
231  ierr = VecDot(arcPtr->dx,arcPtr->dx,&arcPtr->dx2); CHKERRG(ierr);
232  PetscPrintf(PETSC_COMM_WORLD,"\tdlambda = %6.4e dx2 = %6.4e\n",arcPtr->dLambda,arcPtr->dx2);
233  MoFEMFunctionReturnHot(0);
234  }
235 
236  MoFEMErrorCode calculate_init_dlambda(double *dlambda) {
237  MoFEMFunctionBeginHot;
238 
239  *dlambda = arcPtr->s/(arcPtr->beta*sqrt(arcPtr->F_lambda2));
240  PetscPrintf(
241  PETSC_COMM_WORLD,
242  "\tInit dlambda = %6.4e s = %6.4e beta = %6.4e F_lambda2 = %6.4e\n",
243  *dlambda,arcPtr->s,arcPtr->beta,arcPtr->F_lambda2
244  );
245  double a = *dlambda;
246  if(a - a != 0) {
247  std::ostringstream sss;
248  sss << "s " << arcPtr->s << " " << arcPtr->beta << " " << arcPtr->F_lambda2;
249  SETERRQ(PETSC_COMM_SELF,1,sss.str().c_str());
250  }
251 
252  MoFEMFunctionReturnHot(0);
253  }
254 
255  MoFEMErrorCode set_dlambda_to_x(Vec &x,double dlambda) {
256  MoFEMFunctionBeginHot;
257 
258  if(arcPtr->getPetscLocalDofIdx()!=-1) {
259  double *array;
260  ierr = VecGetArray(x,&array); CHKERRG(ierr);
261  double lambda_old = array[arcPtr->getPetscLocalDofIdx()];
262  if(!(dlambda == dlambda)) {
263  std::ostringstream sss;
264  sss << "s " << arcPtr->s << " " << arcPtr->beta << " " << arcPtr->F_lambda2;
265  SETERRQ(PETSC_COMM_SELF,1,sss.str().c_str());
266  }
267  array[arcPtr->getPetscLocalDofIdx()] = lambda_old + dlambda;
268  PetscPrintf(PETSC_COMM_WORLD,"\tlambda = %6.4e, %6.4e (%6.4e)\n",
269  lambda_old, array[arcPtr->getPetscLocalDofIdx()], dlambda);
270  ierr = VecRestoreArray(x,&array); CHKERRG(ierr);
271  }
272 
273  MoFEMFunctionReturnHot(0);
274  }
275 
276 };
277 
278 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
CohesiveElement::ArcLengthIntElemFEMethod::calculate_db
virtual MoFEMErrorCode calculate_db()
Definition: InterfaceGapArcLengthControl.hpp:145
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
CohesiveElement::ArcLengthIntElemFEMethod::Edges4
Range Edges4
Definition: InterfaceGapArcLengthControl.hpp:19
CohesiveElement::ArcLengthIntElemFEMethod
Definition: InterfaceGapArcLengthControl.hpp:12
EntityHandle
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
CohesiveElement::ArcLengthIntElemFEMethod::ArcLengthIntElemFEMethod
ArcLengthIntElemFEMethod(moab::Interface &moab, boost::shared_ptr< ArcLengthCtx > arcptr)
Definition: InterfaceGapArcLengthControl.hpp:24
MOAB_THROW
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:554
CohesiveElement::ArcLengthIntElemFEMethod::operator()
MoFEMErrorCode operator()()
Definition: InterfaceGapArcLengthControl.hpp:172
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
CohesiveElement::ArcLengthIntElemFEMethod::Nodes4
Range Nodes4
Definition: InterfaceGapArcLengthControl.hpp:20
CohesiveElement::ArcLengthIntElemFEMethod::GhostLambdaInt
Vec GhostLambdaInt
Definition: InterfaceGapArcLengthControl.hpp:17
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
CohesiveElement::ArcLengthIntElemFEMethod::calculate_init_dlambda
MoFEMErrorCode calculate_init_dlambda(double *dlambda)
Definition: InterfaceGapArcLengthControl.hpp:236
FEMethod
a
constexpr double a
Definition: approx_sphere.cpp:30
CohesiveElement::ArcLengthIntElemFEMethod::calculate_dx_and_dlambda
MoFEMErrorCode calculate_dx_and_dlambda(Vec &x)
Definition: InterfaceGapArcLengthControl.hpp:214
CohesiveElement::ArcLengthIntElemFEMethod::Nodes3
Range Nodes3
Definition: InterfaceGapArcLengthControl.hpp:20
CohesiveElement::ArcLengthIntElemFEMethod::set_dlambda_to_x
MoFEMErrorCode set_dlambda_to_x(Vec &x, double dlambda)
Definition: InterfaceGapArcLengthControl.hpp:255
CohesiveElement::ArcLengthIntElemFEMethod::~ArcLengthIntElemFEMethod
~ArcLengthIntElemFEMethod()
Definition: InterfaceGapArcLengthControl.hpp:59
CohesiveElement::ArcLengthIntElemFEMethod::remove_damaged_prisms_nodes
MoFEMErrorCode remove_damaged_prisms_nodes()
remove nodes of prims which are fully damaged
Definition: InterfaceGapArcLengthControl.hpp:66
CohesiveElement::ArcLengthIntElemFEMethod::thDamagedPrism
Tag thDamagedPrism
Definition: InterfaceGapArcLengthControl.hpp:22
CohesiveElement::ArcLengthIntElemFEMethod::Faces3
Range Faces3
Definition: InterfaceGapArcLengthControl.hpp:18
Range
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
CohesiveElement::ArcLengthIntElemFEMethod::postProcess
MoFEMErrorCode postProcess()
Definition: InterfaceGapArcLengthControl.hpp:197
CohesiveElement::ArcLengthIntElemFEMethod::lambda_int
double lambda_int
Definition: InterfaceGapArcLengthControl.hpp:87
CohesiveElement
Definition: arc_length_interface.cpp:29
CohesiveElement::ArcLengthIntElemFEMethod::arcPtr
boost::shared_ptr< ArcLengthCtx > arcPtr
Definition: InterfaceGapArcLengthControl.hpp:15
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
CohesiveElement::ArcLengthIntElemFEMethod::Faces4
Range Faces4
Definition: InterfaceGapArcLengthControl.hpp:18
CohesiveElement::ArcLengthIntElemFEMethod::mOab
moab::Interface & mOab
Definition: InterfaceGapArcLengthControl.hpp:13
CohesiveElement::ArcLengthIntElemFEMethod::preProcess
MoFEMErrorCode preProcess()
Definition: InterfaceGapArcLengthControl.hpp:88
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
CohesiveElement::ArcLengthIntElemFEMethod::Edges3
Range Edges3
Definition: InterfaceGapArcLengthControl.hpp:19
CohesiveElement::ArcLengthIntElemFEMethod::calculate_lambda_int
MoFEMErrorCode calculate_lambda_int(double &_lambda_int_)
Definition: InterfaceGapArcLengthControl.hpp:103
Generated by Doxygen 1.8.17 and hosted at University of Glasgow