v0.15.0
Loading...
Searching...
No Matches
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
10namespace CohesiveElement {
11
12struct 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}
a
constexpr double a
Definition approx_sphere.cpp:30
ierr
static PetscErrorCode ierr
Definition base_functions.c:13
MOAB_THROW
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition definitions.h:557
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition definitions.h:228
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition definitions.h:499
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
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::calculate_lambda_int
MoFEMErrorCode calculate_lambda_int(double &_lambda_int_)
Definition InterfaceGapArcLengthControl.hpp:103
CohesiveElement::ArcLengthIntElemFEMethod::ArcLengthIntElemFEMethod
ArcLengthIntElemFEMethod(moab::Interface &moab, boost::shared_ptr< ArcLengthCtx > arcptr)
Definition InterfaceGapArcLengthControl.hpp:24
CohesiveElement::ArcLengthIntElemFEMethod::set_dlambda_to_x
MoFEMErrorCode set_dlambda_to_x(Vec &x, double dlambda)
Definition InterfaceGapArcLengthControl.hpp:255
Generated by Doxygen 1.12.0 and hosted at University of Glasgow