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