v0.13.2
Loading...
Searching...
No Matches
EshelbianPlasticity.cpp
Go to the documentation of this file.
1/**
2 * \file EshelbianPlasticity.cpp
3 * \example EshelbianPlasticity.cpp
4 *
5 * \brief Eshelbian plasticity implementation
6 */
7
8#include <MoFEM.hpp>
9using namespace MoFEM;
10
11#include <BasicFiniteElements.hpp>
12#include <CGGTonsorialBubbleBase.hpp>
13
14#include <EshelbianPlasticity.hpp>
15#include <boost/math/constants/constants.hpp>
16
17namespace EshelbianPlasticity {
18
19MoFEMErrorCode
20EshelbianCore::query_interface(boost::typeindex::type_index type_index,
21 UnknownInterface **iface) const {
22 *iface = const_cast<EshelbianCore *>(this);
23 return 0;
24}
25
26MoFEMErrorCode OpJacobian::doWork(int side, EntityType type, EntData &data) {
27 MoFEMFunctionBegin;
28 if (type != MBVERTEX)
29 MoFEMFunctionReturnHot(0);
30
31 if (evalRhs)
32 CHKERR evaluateRhs(data);
33
34 if (evalLhs)
35 CHKERR evaluateLhs(data);
36
37 MoFEMFunctionReturn(0);
38}
39
40EshelbianCore::EshelbianCore(MoFEM::Interface &m_field)
41 : mField(m_field), piolaStress("P"), eshelbyStress("S"), spatialDisp("w"),
42 materialDisp("W"), streachTensor("u"), rotAxis("omega"),
43 materialGradient("G"), tauField("TAU"), lambdaField("LAMBDA"),
44 bubbleField("BUBBLE"), elementVolumeName("EP"),
45 naturalBcElement("NATURAL_BC"), essentialBcElement("ESSENTIAL_BC") {
46
47 ierr = getOptions();
48 CHKERRABORT(PETSC_COMM_WORLD, ierr);
49}
50
51EshelbianCore::~EshelbianCore() {}
52
53MoFEMErrorCode EshelbianCore::getOptions() {
54 MoFEMFunctionBegin;
55 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Eshelbian plasticity",
56 "none");
57
58 spaceOrder = 1;
59 CHKERR PetscOptionsInt("-space_order", "approximation oder for space", "",
60 spaceOrder, &spaceOrder, PETSC_NULL);
61 materialOrder = 1;
62 CHKERR PetscOptionsInt("-material_order", "approximation oder for material",
63 "", materialOrder, &materialOrder, PETSC_NULL);
64
65 alphaU = 0;
66 CHKERR PetscOptionsScalar("-viscosity_alpha_u", "viscosity", "", alphaU,
67 &alphaU, PETSC_NULL);
68
69 alphaW = 0;
70 CHKERR PetscOptionsScalar("-viscosity_alpha_w", "viscosity", "", alphaW,
71 &alphaW, PETSC_NULL);
72
73 alphaRho = 0;
74 CHKERR PetscOptionsScalar("-density_alpha_rho", "density", "", alphaRho,
75 &alphaRho, PETSC_NULL);
76
77 precEps = 1e-3;
78 CHKERR PetscOptionsScalar("-preconditioner_eps", "preconditioner_eps", "",
79 precEps, &precEps, PETSC_NULL);
80
81 ierr = PetscOptionsEnd();
82 CHKERRG(ierr);
83 MoFEMFunctionReturn(0);
84}
85
86MoFEMErrorCode EshelbianCore::addFields(const EntityHandle meshset) {
87 MoFEMFunctionBegin;
88
89 Range tets;
90 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
91 Range tets_skin_part;
92 Skinner skin(&mField.get_moab());
93 CHKERR skin.find_skin(0, tets, false, tets_skin_part);
94 ParallelComm *pcomm =
95 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
96 Range tets_skin;
97 CHKERR pcomm->filter_pstatus(tets_skin_part,
98 PSTATUS_SHARED | PSTATUS_MULTISHARED,
99 PSTATUS_NOT, -1, &tets_skin);
100
101 auto subtract_faces_where_displacements_are_applied =
102 [&](const std::string disp_block_set_name) {
103 MoFEMFunctionBegin;
104 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
105 if (it->getName().compare(0, disp_block_set_name.length(),
106 disp_block_set_name) == 0) {
107 Range faces;
108 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2,
109 faces, true);
110 tets_skin = subtract(tets_skin, faces);
111 }
112 }
113 MoFEMFunctionReturn(0);
114 };
115 CHKERR subtract_faces_where_displacements_are_applied("SPATIAL_DISP_BC");
116 CHKERR subtract_faces_where_displacements_are_applied("SPATIAL_ROTATION_BC");
117 CHKERR subtract_faces_where_displacements_are_applied("SPATIAL_TRACTION_BC");
118
119 Range faces;
120 CHKERR mField.get_moab().get_adjacencies(tets, 2, true, faces,
121 moab::Interface::UNION);
122 Range faces_not_on_the_skin = subtract(faces, tets_skin);
123
124 auto add_hdiv_field = [&](const std::string field_name, const int order,
125 const int dim) {
126 MoFEMFunctionBegin;
127 CHKERR mField.add_field(field_name, HDIV, DEMKOWICZ_JACOBI_BASE, dim,
128 MB_TAG_SPARSE, MF_ZERO);
129 CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
130 CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
131 CHKERR mField.set_field_order(faces_not_on_the_skin, field_name, order);
132 CHKERR mField.set_field_order(tets_skin, field_name, 0);
133 MoFEMFunctionReturn(0);
134 };
135
136 auto add_hdiv_rt_field = [&](const std::string field_name, const int order,
137 const int dim) {
138 MoFEMFunctionBegin;
139 CHKERR mField.add_field(field_name, HDIV, DEMKOWICZ_JACOBI_BASE, dim,
140 MB_TAG_DENSE, MF_ZERO);
141 CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
142 CHKERR mField.set_field_order(meshset, MBTET, field_name, 0);
143 CHKERR mField.set_field_order(tets_skin, field_name, order);
144 MoFEMFunctionReturn(0);
145 };
146
147 auto add_l2_field = [this, meshset](const std::string field_name,
148 const int order, const int dim) {
149 MoFEMFunctionBegin;
150 CHKERR mField.add_field(field_name, L2, AINSWORTH_LEGENDRE_BASE, dim,
151 MB_TAG_DENSE, MF_ZERO);
152 CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
153 CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
154 MoFEMFunctionReturn(0);
155 };
156
157 auto add_h1_field = [this, meshset](const std::string field_name,
158 const int order, const int dim) {
159 MoFEMFunctionBegin;
160 CHKERR mField.add_field(field_name, H1, AINSWORTH_LEGENDRE_BASE, dim,
161 MB_TAG_DENSE, MF_ZERO);
162 CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
163 CHKERR mField.set_field_order(meshset, MBVERTEX, field_name, 1);
164 CHKERR mField.set_field_order(meshset, MBEDGE, field_name, order);
165 CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
166 CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
167 MoFEMFunctionReturn(0);
168 };
169
170 auto add_bubble_field = [this, meshset](const std::string field_name,
171 const int order, const int dim) {
172 MoFEMFunctionBegin;
173 CHKERR mField.add_field(field_name, HDIV, USER_BASE, dim, MB_TAG_DENSE,
174 MF_ZERO);
175 // Modify field
176 auto field_ptr = mField.get_field_structure(field_name);
177 auto field_order_table =
178 const_cast<Field *>(field_ptr)->getFieldOrderTable();
179 auto get_cgg_bubble_order_zero = [](int p) { return 0; };
180 auto get_cgg_bubble_order_tet = [](int p) {
181 return NBVOLUMETET_CCG_BUBBLE(p);
182 };
183 field_order_table[MBVERTEX] = get_cgg_bubble_order_zero;
184 field_order_table[MBEDGE] = get_cgg_bubble_order_zero;
185 field_order_table[MBTRI] = get_cgg_bubble_order_zero;
186 field_order_table[MBTET] = get_cgg_bubble_order_tet;
187 CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
188 CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
189 CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
190 MoFEMFunctionReturn(0);
191 };
192
193 // spatial fields
194 CHKERR add_hdiv_field(piolaStress, spaceOrder, 3);
195 CHKERR add_bubble_field(bubbleField, spaceOrder, 1);
196 CHKERR add_l2_field(spatialDisp, spaceOrder - 1, 3);
197 CHKERR add_l2_field(rotAxis, spaceOrder - 1, 3);
198 CHKERR add_l2_field(streachTensor, spaceOrder, 6);
199
200 // material fields
201 CHKERR add_hdiv_field(eshelbyStress, materialOrder, 3);
202 CHKERR add_l2_field(materialGradient, materialOrder - 1, 9);
203 // CHKERR add_l2_field(materialDisp, materialOrder - 1, 3);
204 // CHKERR add_l2_field(tauField, materialOrder - 1, 1);
205 // CHKERR add_l2_field(lambdaField, materialOrder - 1, 1);
206
207 // Add history filedes
208 CHKERR add_l2_field(materialGradient + "0", materialOrder - 1, 9);
209
210 CHKERR mField.build_fields();
211
212 MoFEMFunctionReturn(0);
213}
214
215MoFEMErrorCode
216EshelbianCore::addVolumeFiniteElement(const EntityHandle meshset) {
217 MoFEMFunctionBegin;
218
219 // set finite element fields
220 auto add_field_to_fe = [this](const std::string fe,
221 const std::string field_name) {
222 MoFEMFunctionBegin;
223 CHKERR mField.modify_finite_element_add_field_row(fe, field_name);
224 CHKERR mField.modify_finite_element_add_field_col(fe, field_name);
225 CHKERR mField.modify_finite_element_add_field_data(fe, field_name);
226 MoFEMFunctionReturn(0);
227 };
228
229 if (!mField.check_finite_element(elementVolumeName)) {
230 CHKERR mField.add_finite_element(elementVolumeName, MF_ZERO);
231 CHKERR mField.add_ents_to_finite_element_by_type(meshset, MBTET,
232 elementVolumeName);
233
234 CHKERR add_field_to_fe(elementVolumeName, piolaStress);
235 CHKERR add_field_to_fe(elementVolumeName, bubbleField);
236 CHKERR add_field_to_fe(elementVolumeName, eshelbyStress);
237 CHKERR add_field_to_fe(elementVolumeName, streachTensor);
238 CHKERR add_field_to_fe(elementVolumeName, rotAxis);
239 CHKERR add_field_to_fe(elementVolumeName, spatialDisp);
240 CHKERR add_field_to_fe(elementVolumeName, streachTensor);
241 CHKERR add_field_to_fe(elementVolumeName, materialGradient);
242
243 CHKERR mField.modify_finite_element_add_field_data(elementVolumeName,
244 materialGradient + "0");
245 }
246
247 // build finite elements data structures
248 CHKERR mField.build_finite_elements(elementVolumeName);
249
250 MoFEMFunctionReturn(0);
251}
252
253MoFEMErrorCode
254EshelbianCore::addBoundaryFiniteElement(const EntityHandle meshset) {
255 MoFEMFunctionBegin;
256
257 auto bc_elements_add_to_range = [&](const std::string disp_block_set_name,
258 Range &r) {
259 MoFEMFunctionBegin;
260 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
261 if (it->getName().compare(0, disp_block_set_name.length(),
262 disp_block_set_name) == 0) {
263 Range faces;
264 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
265 true);
266 r.merge(faces);
267 }
268 }
269 MoFEMFunctionReturn(0);
270 };
271
272 // set finite element fields
273 auto add_field_to_fe = [this](const std::string fe,
274 const std::string field_name) {
275 MoFEMFunctionBegin;
276 CHKERR mField.modify_finite_element_add_field_row(fe, field_name);
277 CHKERR mField.modify_finite_element_add_field_col(fe, field_name);
278 CHKERR mField.modify_finite_element_add_field_data(fe, field_name);
279 MoFEMFunctionReturn(0);
280 };
281
282 Range natural_bc_elements;
283 CHKERR bc_elements_add_to_range("SPATIAL_DISP_BC", natural_bc_elements);
284 CHKERR bc_elements_add_to_range("SPATIAL_ROTATION_BC", natural_bc_elements);
285 Range essentail_bc_elements;
286 CHKERR bc_elements_add_to_range("SPATIAL_TRACTION_BC", essentail_bc_elements);
287
288 CHKERR mField.add_finite_element(naturalBcElement, MF_ZERO);
289 CHKERR mField.add_ents_to_finite_element_by_type(natural_bc_elements, MBTRI,
290 naturalBcElement);
291 CHKERR add_field_to_fe(naturalBcElement, piolaStress);
292 CHKERR add_field_to_fe(naturalBcElement, eshelbyStress);
293 CHKERR mField.build_finite_elements(naturalBcElement);
294
295 CHKERR mField.add_finite_element(essentialBcElement, MF_ZERO);
296 CHKERR mField.add_ents_to_finite_element_by_type(essentail_bc_elements, MBTRI,
297 essentialBcElement);
298 CHKERR add_field_to_fe(essentialBcElement, piolaStress);
299 CHKERR add_field_to_fe(essentialBcElement, eshelbyStress);
300 CHKERR mField.build_finite_elements(essentialBcElement);
301
302 MoFEMFunctionReturn(0);
303}
304
305MoFEMErrorCode EshelbianCore::addDMs(const BitRefLevel bit) {
306 MoFEMFunctionBegin;
307
308 // find adjacencies between finite elements and dofs
309 CHKERR mField.build_adjacencies(bit, QUIET);
310
311 // Create coupled problem
312 dM = createSmartDM(mField.get_comm(), "DMMOFEM");
313 CHKERR DMMoFEMCreateMoFEM(dM, &mField, "ESHELBY_PLASTICITY", bit,
314 BitRefLevel().set());
315 CHKERR DMMoFEMSetDestroyProblem(dM, PETSC_TRUE);
316 CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
317 CHKERR DMMoFEMAddElement(dM, elementVolumeName.c_str());
318 CHKERR DMMoFEMAddElement(dM, naturalBcElement.c_str());
319 CHKERR DMMoFEMAddElement(dM, essentialBcElement.c_str());
320 mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
321 CHKERR DMSetUp(dM);
322 mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
323
324 auto remove_dofs_on_essential_spatial_stress_boundary =
325 [&](const std::string prb_name) {
326 MoFEMFunctionBegin;
327 for (int d : {0, 1, 2})
328 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
329 prb_name, piolaStress, (*bcSpatialFreeTraction)[d], d, d, 0, 100,
330 NOISY, true);
331 MoFEMFunctionReturn(0);
332 };
333 CHKERR remove_dofs_on_essential_spatial_stress_boundary("ESHELBY_PLASTICITY");
334
335 // Create elastic sub-problem
336 dmElastic = createSmartDM(mField.get_comm(), "DMMOFEM");
337 CHKERR DMMoFEMCreateSubDM(dmElastic, dM, "ELASTIC_PROBLEM");
338 CHKERR DMMoFEMSetDestroyProblem(dmElastic, PETSC_TRUE);
339 CHKERR DMMoFEMAddSubFieldRow(dmElastic, piolaStress.c_str());
340 CHKERR DMMoFEMAddSubFieldRow(dmElastic, bubbleField.c_str());
341 CHKERR DMMoFEMAddSubFieldRow(dmElastic, streachTensor.c_str());
342 CHKERR DMMoFEMAddSubFieldRow(dmElastic, rotAxis.c_str());
343 CHKERR DMMoFEMAddSubFieldRow(dmElastic, spatialDisp.c_str());
344 CHKERR DMMoFEMAddElement(dmElastic, elementVolumeName.c_str());
345 CHKERR DMMoFEMAddElement(dmElastic, naturalBcElement.c_str());
346 CHKERR DMMoFEMAddElement(dmElastic, essentialBcElement.c_str());
347 CHKERR DMMoFEMSetSquareProblem(dmElastic, PETSC_TRUE);
348 CHKERR DMMoFEMSetIsPartitioned(dmElastic, PETSC_TRUE);
349 CHKERR DMSetUp(dmElastic);
350
351 // Create elastic streach-problem
352 dmElasticSchurStreach = createSmartDM(mField.get_comm(), "DMMOFEM");
353 CHKERR DMMoFEMCreateSubDM(dmElasticSchurStreach, dmElastic,
354 "ELASTIC_PROBLEM_STREACH_SCHUR");
355 CHKERR DMMoFEMSetDestroyProblem(dmElasticSchurStreach, PETSC_TRUE);
356 CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurStreach, piolaStress.c_str());
357 CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurStreach, bubbleField.c_str());
358 CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurStreach, rotAxis.c_str());
359 CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurStreach, spatialDisp.c_str());
360 CHKERR DMMoFEMAddElement(dmElasticSchurStreach, elementVolumeName.c_str());
361 CHKERR DMMoFEMAddElement(dmElasticSchurStreach, naturalBcElement.c_str());
362 CHKERR DMMoFEMAddElement(dmElasticSchurStreach, essentialBcElement.c_str());
363 CHKERR DMMoFEMSetSquareProblem(dmElasticSchurStreach, PETSC_TRUE);
364 CHKERR DMMoFEMSetIsPartitioned(dmElasticSchurStreach, PETSC_TRUE);
365 CHKERR DMSetUp(dmElasticSchurStreach);
366
367 // Create elastic bubble-problem
368 dmElasticSchurBubble = createSmartDM(mField.get_comm(), "DMMOFEM");
369 CHKERR DMMoFEMCreateSubDM(dmElasticSchurBubble, dmElasticSchurStreach,
370 "ELASTIC_PROBLEM_BUBBLE_SCHUR");
371 CHKERR DMMoFEMSetDestroyProblem(dmElasticSchurBubble, PETSC_TRUE);
372 CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurBubble, piolaStress.c_str());
373 CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurBubble, rotAxis.c_str());
374 CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurBubble, spatialDisp.c_str());
375 CHKERR DMMoFEMAddElement(dmElasticSchurBubble, elementVolumeName.c_str());
376 CHKERR DMMoFEMAddElement(dmElasticSchurBubble, naturalBcElement.c_str());
377 CHKERR DMMoFEMAddElement(dmElasticSchurBubble, essentialBcElement.c_str());
378 CHKERR DMMoFEMSetSquareProblem(dmElasticSchurBubble, PETSC_TRUE);
379 CHKERR DMMoFEMSetIsPartitioned(dmElasticSchurBubble, PETSC_TRUE);
380 CHKERR DMSetUp(dmElasticSchurBubble);
381
382 // Create elastic omega-problem
383 dmElasticSchurOmega = createSmartDM(mField.get_comm(), "DMMOFEM");
384 CHKERR DMMoFEMCreateSubDM(dmElasticSchurOmega, dmElasticSchurBubble,
385 "ELASTIC_PROBLEM_OMEGA_SCHUR");
386 CHKERR DMMoFEMSetDestroyProblem(dmElasticSchurOmega, PETSC_TRUE);
387 CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurOmega, piolaStress.c_str());
388 CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurOmega, spatialDisp.c_str());
389 CHKERR DMMoFEMAddElement(dmElasticSchurOmega, elementVolumeName.c_str());
390 CHKERR DMMoFEMAddElement(dmElasticSchurOmega, naturalBcElement.c_str());
391 CHKERR DMMoFEMAddElement(dmElasticSchurOmega, essentialBcElement.c_str());
392 CHKERR DMMoFEMSetSquareProblem(dmElasticSchurOmega, PETSC_TRUE);
393 CHKERR DMMoFEMSetIsPartitioned(dmElasticSchurOmega, PETSC_TRUE);
394 CHKERR DMSetUp(dmElasticSchurOmega);
395
396 // Create elastic tet_stress-problem
397 dmElasticSchurSpatialDisp = createSmartDM(mField.get_comm(), "DMMOFEM");
398 CHKERR DMMoFEMCreateSubDM(dmElasticSchurSpatialDisp, dmElasticSchurOmega,
399 "ELASTIC_PROBLEM_SPATIAL_DISP_SCHUR");
400 CHKERR DMMoFEMSetDestroyProblem(dmElasticSchurSpatialDisp, PETSC_TRUE);
401 CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurSpatialDisp, piolaStress.c_str());
402 CHKERR DMMoFEMAddElement(dmElasticSchurSpatialDisp,
403 elementVolumeName.c_str());
404 CHKERR DMMoFEMAddElement(dmElasticSchurSpatialDisp, naturalBcElement.c_str());
405 CHKERR DMMoFEMAddElement(dmElasticSchurSpatialDisp,
406 essentialBcElement.c_str());
407 CHKERR DMMoFEMSetSquareProblem(dmElasticSchurSpatialDisp, PETSC_TRUE);
408 CHKERR DMMoFEMSetIsPartitioned(dmElasticSchurSpatialDisp, PETSC_TRUE);
409 CHKERR DMSetUp(dmElasticSchurSpatialDisp);
410
411 {
412 PetscSection section;
413 CHKERR mField.getInterface<ISManager>()->sectionCreate("ELASTIC_PROBLEM",
414 &section);
415 CHKERR DMSetDefaultSection(dmElastic, section);
416 CHKERR DMSetDefaultGlobalSection(dmElastic, section);
417 CHKERR PetscSectionDestroy(&section);
418 }
419
420 MoFEMFunctionReturn(0);
421}
422
423BcDisp::BcDisp(std::string name, std::vector<double> &attr, Range &faces)
424 : blockName(name), faces(faces) {
425 vals.resize(3, false);
426 flags.resize(3, false);
427 for (int ii = 0; ii != 3; ++ii) {
428 vals[ii] = attr[ii];
429 flags[ii] = static_cast<int>(attr[ii + 3]);
430 }
431}
432
433BcRot::BcRot(std::string name, std::vector<double> &attr, Range &faces)
434 : blockName(name), faces(faces) {
435 vals.resize(3, false);
436 for (int ii = 0; ii != 3; ++ii) {
437 vals[ii] = attr[ii];
438 }
439 theta = attr[3];
440}
441
442TractionBc::TractionBc(std::string name, std::vector<double> &attr,
443 Range &faces)
444 : blockName(name), faces(faces) {
445 vals.resize(3, false);
446 flags.resize(3, false);
447 for (int ii = 0; ii != 3; ++ii) {
448 vals[ii] = attr[ii];
449 flags[ii] = static_cast<int>(attr[ii + 3]);
450 }
451}
452
453MoFEMErrorCode
454EshelbianCore::getTractionFreeBc(const EntityHandle meshset,
455 boost::shared_ptr<TractionFreeBc> &bc_ptr,
456 const std::string disp_block_set_name,
457 const std::string rot_block_set_name) {
458 MoFEMFunctionBegin;
459 Range tets;
460 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
461 Range tets_skin_part;
462 Skinner skin(&mField.get_moab());
463 CHKERR skin.find_skin(0, tets, false, tets_skin_part);
464 ParallelComm *pcomm =
465 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
466 Range tets_skin;
467 CHKERR pcomm->filter_pstatus(tets_skin_part,
468 PSTATUS_SHARED | PSTATUS_MULTISHARED,
469 PSTATUS_NOT, -1, &tets_skin);
470
471 bc_ptr->resize(3);
472 for (int dd = 0; dd != 3; ++dd)
473 (*bc_ptr)[dd] = tets_skin;
474
475 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
476 if (it->getName().compare(0, disp_block_set_name.length(),
477 disp_block_set_name) == 0) {
478 std::vector<double> block_attributes;
479 CHKERR it->getAttributes(block_attributes);
480 if (block_attributes.size() != 6) {
481 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
482 "In block %s six attributes are required for given BC "
483 "blockset (3 values + "
484 "3 flags) != %d",
485 it->getName().c_str(), block_attributes.size());
486 }
487 Range faces;
488 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
489 true);
490 if (block_attributes[3] != 0)
491 (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
492 if (block_attributes[4] != 0)
493 (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
494 if (block_attributes[5] != 0)
495 (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
496 }
497 if (it->getName().compare(0, rot_block_set_name.length(),
498 rot_block_set_name) == 0) {
499 Range faces;
500 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
501 true);
502 (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
503 (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
504 (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
505 }
506 }
507
508 // for (int dd = 0; dd != 3; ++dd) {
509 // EntityHandle meshset;
510 // CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
511 // CHKERR mField.get_moab().add_entities(meshset, (*bc_ptr)[dd]);
512 // std::string file_name = disp_block_set_name +
513 // "_traction_free_bc_" + boost::lexical_cast<std::string>(dd) + ".vtk";
514 // CHKERR mField.get_moab().write_file(file_name.c_str(), " VTK ", "",
515 // &meshset, 1);
516 // CHKERR mField.get_moab().delete_entities(&meshset, 1);
517 // }
518
519 MoFEMFunctionReturn(0);
520}
521
522/**
523 * @brief Set integration rule on element
524 * \param order on row
525 * \param order on column
526 * \param order on data
527 *
528 * Use maximal oder on data in order to determine integration rule
529 *
530 */
531struct VolRule {
532 int operator()(int p_row, int p_col, int p_data) const {
533 return 2 * (p_data + 1);
534 }
535};
536
537struct FaceRule {
538 int operator()(int p_row, int p_col, int p_data) const {
539 return 2 * (p_data);
540 }
541};
542
543struct CGGUserPolynomialBase : public BaseFunction {
544
545 CGGUserPolynomialBase() {}
546 ~CGGUserPolynomialBase() {}
547
548 MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
549 BaseFunctionUnknownInterface **iface) const {
550 *iface = const_cast<CGGUserPolynomialBase *>(this);
551 return 0;
552 }
553
554 MoFEMErrorCode getValue(MatrixDouble &pts,
555 boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
556 MoFEMFunctionBeginHot;
557
558 cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
559
560 int nb_gauss_pts = pts.size2();
561 if (!nb_gauss_pts) {
562 MoFEMFunctionReturnHot(0);
563 }
564
565 if (pts.size1() < 3) {
566 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
567 "Wrong dimension of pts, should be at least 3 rows with "
568 "coordinates");
569 }
570
571 switch (cTx->sPace) {
572 case HDIV:
573 CHKERR getValueHdivForCGGBubble(pts);
574 break;
575 default:
576 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
577 }
578
579 MoFEMFunctionReturnHot(0);
580 }
581
582private:
583 EntPolynomialBaseCtx *cTx;
584
585 MatrixDouble shapeFun;
586
587 MoFEMErrorCode getValueHdivForCGGBubble(MatrixDouble &pts) {
588 MoFEMFunctionBegin;
589
590 // This should be used only in case USER_BASE is selected
591 if (cTx->bAse != USER_BASE) {
592 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
593 "Wrong base, should be USER_BASE");
594 }
595 // get access to data structures on element
596 DataForcesAndSourcesCore &data = cTx->dAta;
597 // Get approximation order on element. Note that bubble functions are only
598 // on tetrahedron.
599 const int order = data.dataOnEntities[MBTET][0].getOrder();
600 /// number of integration points
601 const int nb_gauss_pts = pts.size2();
602 // number of base functions
603
604 // calculate shape functions, i.e. barycentric coordinates
605 shapeFun.resize(nb_gauss_pts, 4, false);
606 CHKERR ShapeMBTET(&*shapeFun.data().begin(), &pts(0, 0), &pts(1, 0),
607 &pts(2, 0), nb_gauss_pts);
608 // direvatives of shape functions
609 double diff_shape_fun[12];
610 CHKERR ShapeDiffMBTET(diff_shape_fun);
611
612 const int nb_base_functions = NBVOLUMETET_CCG_BUBBLE(order);
613 // get base functions and set size
614 MatrixDouble &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
615 phi.resize(nb_gauss_pts, 9 * nb_base_functions, false);
616 // finally calculate base functions
617 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_phi(
618 &phi(0, 0), &phi(0, 1), &phi(0, 2),
619
620 &phi(0, 3), &phi(0, 4), &phi(0, 5),
621
622 &phi(0, 6), &phi(0, 7), &phi(0, 8));
623 CHKERR CGG_BubbleBase_MBTET(order, &shapeFun(0, 0), diff_shape_fun, t_phi,
624 nb_gauss_pts);
625
626 MoFEMFunctionReturn(0);
627 }
628};
629
630MoFEMErrorCode EshelbianCore::setBaseVolumeElementOps(
631 const int tag, const bool do_rhs, const bool do_lhs,
632 boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe) {
633 MoFEMFunctionBegin;
634 fe = boost::make_shared<EpElement<VolumeElementForcesAndSourcesCore>>(mField);
635
636 fe->getUserPolynomialBase() =
637 boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
638 fe->getOpPtrVector().push_back(new OpL2Transform());
639
640 // set integration rule
641 fe->getRuleHook = VolRule();
642
643 if (!dataAtPts) {
644 dataAtPts =
645 boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
646 dataAtPts->physicsPtr = physicalEquations;
647 }
648
649 // calculate fields values
650 fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
651 piolaStress, dataAtPts->getApproxPAtPts()));
652 fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
653 bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
654 fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
655 piolaStress, dataAtPts->getDivPAtPts()));
656 fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
657 eshelbyStress, dataAtPts->getApproxSigmaAtPts()));
658 fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
659 eshelbyStress, dataAtPts->getDivSigmaAtPts()));
660 fe->getOpPtrVector().push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
661 streachTensor, dataAtPts->getLogStreachTensorAtPts(), MBTET));
662 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
663 rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
664 fe->getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
665 materialGradient, dataAtPts->getBigGAtPts(), MBTET));
666 fe->getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
667 materialGradient + "0", dataAtPts->getBigG0AtPts(), MBTET));
668 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
669 spatialDisp, dataAtPts->getSmallWAtPts(), MBTET));
670
671 // velocities
672 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
673 spatialDisp, dataAtPts->getSmallWDotAtPts(), MBTET));
674 fe->getOpPtrVector().push_back(
675 new OpCalculateTensor2SymmetricFieldValuesDot<3>(
676 streachTensor, dataAtPts->getLogStreachDotTensorAtPts(), MBTET));
677 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
678 rotAxis, dataAtPts->getRotAxisDotAtPts(), MBTET));
679
680 // acceleration
681 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
682 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDotDot<3>(
683 spatialDisp, dataAtPts->getSmallWDotDotAtPts(), MBTET));
684 }
685
686 // calculate other derived quantities
687 fe->getOpPtrVector().push_back(
688 new OpCalculateRotationAndSpatialGradient(rotAxis, dataAtPts));
689
690 // evaluate integration points
691 fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
692 spatialDisp, tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
693
694 MoFEMFunctionReturn(0);
695}
696
697MoFEMErrorCode EshelbianCore::setGenericVolumeElementOps(
698 const int tag, const bool add_elastic, const bool add_material,
699 boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe_rhs,
700 boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe_lhs) {
701 MoFEMFunctionBegin;
702
703 // Right hand side
704 CHKERR setBaseVolumeElementOps(tag, true, false, fe_rhs);
705
706 // elastic
707 if (add_elastic) {
708 fe_rhs->getOpPtrVector().push_back(
709 new OpSpatialEquilibrium(spatialDisp, dataAtPts, alphaW, alphaRho));
710 fe_rhs->getOpPtrVector().push_back(
711 new OpSpatialRotation(rotAxis, dataAtPts));
712 fe_rhs->getOpPtrVector().push_back(
713 new OpSpatialPhysical(streachTensor, dataAtPts, alphaU));
714 fe_rhs->getOpPtrVector().push_back(
715 new OpSpatialConsistencyP(piolaStress, dataAtPts));
716 fe_rhs->getOpPtrVector().push_back(
717 new OpSpatialConsistencyBubble(bubbleField, dataAtPts));
718 fe_rhs->getOpPtrVector().push_back(
719 new OpSpatialConsistencyDivTerm(piolaStress, dataAtPts));
720 }
721
722 // Left hand side
723 CHKERR setBaseVolumeElementOps(tag, true, true, fe_lhs);
724
725 // elastic
726 if (add_elastic) {
727
728 // Schur
729 fe_lhs->getOpPtrVector().push_back(
730 new OpSpatialSchurBegin(spatialDisp, dataAtPts));
731
732 fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_du(
733 streachTensor, streachTensor, dataAtPts, alphaU));
734 fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dBubble(
735 streachTensor, bubbleField, dataAtPts, true));
736
737 fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dP(
738 spatialDisp, piolaStress, dataAtPts, true));
739 fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dw(
740 spatialDisp, spatialDisp, dataAtPts, alphaW, alphaRho));
741
742 fe_lhs->getOpPtrVector().push_back(
743 new OpSpatialConsistency_dP_domega(piolaStress, rotAxis, dataAtPts));
744 fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
745 bubbleField, rotAxis, dataAtPts));
746
747 fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dP(
748 streachTensor, piolaStress, dataAtPts, true));
749 fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
750 streachTensor, rotAxis, dataAtPts, false));
751
752 fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dP(
753 rotAxis, piolaStress, dataAtPts, false));
754 fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dBubble(
755 rotAxis, bubbleField, dataAtPts, false));
756 fe_lhs->getOpPtrVector().push_back(
757 new OpSpatialRotation_domega_domega(rotAxis, rotAxis, dataAtPts));
758
759 // Schur
760 dataAtPts->ooMatPtr = boost::make_shared<MatrixDouble>();
761 fe_lhs->getOpPtrVector().push_back(
762 new OpSpatialPreconditionMass(rotAxis, dataAtPts->ooMatPtr));
763 if (alphaW < std::numeric_limits<double>::epsilon() &&
764 alphaRho < std::numeric_limits<double>::epsilon()) {
765 dataAtPts->wwMatPtr = boost::make_shared<MatrixDouble>();
766 fe_lhs->getOpPtrVector().push_back(
767 new OpSpatialPreconditionMass(spatialDisp, dataAtPts->wwMatPtr));
768 } else {
769 dataAtPts->wwMatPtr.reset();
770 }
771 fe_lhs->getOpPtrVector().push_back(
772 new OpSpatialSchurEnd(spatialDisp, dataAtPts, precEps));
773
774 if (add_material) {
775 }
776 }
777
778 MoFEMFunctionReturn(0);
779}
780
781MoFEMErrorCode EshelbianCore::setGenericFaceElementOps(
782 const bool add_elastic, const bool add_material,
783 boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> &fe_rhs,
784 boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> &fe_lhs) {
785 MoFEMFunctionBegin;
786
787 fe_rhs =
788 boost::make_shared<EpElement<FaceElementForcesAndSourcesCore>>(mField);
789 fe_lhs =
790 boost::make_shared<EpElement<FaceElementForcesAndSourcesCore>>(mField);
791
792 // set integration rule
793 fe_rhs->getRuleHook = FaceRule();
794 fe_lhs->getRuleHook = FaceRule();
795
796 fe_rhs->getOpPtrVector().push_back(
797 new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
798 fe_lhs->getOpPtrVector().push_back(
799 new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
800
801 if (add_elastic) {
802 fe_rhs->getOpPtrVector().push_back(
803 new OpDispBc(piolaStress, dataAtPts, bcSpatialDispVecPtr));
804 fe_rhs->getOpPtrVector().push_back(
805 new OpRotationBc(piolaStress, dataAtPts, bcSpatialRotationVecPtr));
806 }
807
808 MoFEMFunctionReturn(0);
809}
810
811MoFEMErrorCode EshelbianCore::setElasticElementOps(const int tag) {
812 MoFEMFunctionBegin;
813 CHKERR setGenericVolumeElementOps(tag, true, false, elasticFeRhs,
814 elasticFeLhs);
815 CHKERR setGenericFaceElementOps(true, false, elasticBcRhs, elasticBcLhs);
816 MoFEMFunctionReturn(0);
817}
818
819MoFEMErrorCode EshelbianCore::setElasticElementToTs(DM dm) {
820 MoFEMFunctionBegin;
821 boost::shared_ptr<FEMethod> null;
822 boost::shared_ptr<EpElement<FeTractionBc>> spatial_traction_bc(
823 new EpElement<FeTractionBc>(mField, piolaStress, bcSpatialTraction));
824 schurAssembly = boost::make_shared<EpFEMethod>();
825
826 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
827 CHKERR DMMoFEMTSSetI2Function(dm, DM_NO_ELEMENT, null, spatial_traction_bc,
828 null);
829 CHKERR DMMoFEMTSSetI2Function(dm, elementVolumeName, elasticFeRhs, null,
830 null);
831 CHKERR DMMoFEMTSSetI2Function(dm, naturalBcElement, elasticBcRhs, null,
832 null);
833 CHKERR DMMoFEMTSSetI2Function(dm, DM_NO_ELEMENT, null, null,
834 spatial_traction_bc);
835 CHKERR DMMoFEMTSSetI2Jacobian(dm, DM_NO_ELEMENT, null, schurAssembly, null);
836 CHKERR DMMoFEMTSSetI2Jacobian(dm, elementVolumeName, elasticFeLhs, null,
837 null);
838 CHKERR DMMoFEMTSSetI2Jacobian(dm, naturalBcElement, elasticBcLhs, null,
839 null);
840 CHKERR DMMoFEMTSSetI2Jacobian(dm, DM_NO_ELEMENT, null, null, schurAssembly);
841 } else {
842 CHKERR DMMoFEMTSSetIFunction(dm, DM_NO_ELEMENT, null, spatial_traction_bc,
843 null);
844 CHKERR DMMoFEMTSSetIFunction(dm, elementVolumeName, elasticFeRhs, null,
845 null);
846 CHKERR DMMoFEMTSSetIFunction(dm, naturalBcElement, elasticBcRhs, null,
847 null);
848 CHKERR DMMoFEMTSSetIFunction(dm, DM_NO_ELEMENT, null, null,
849 spatial_traction_bc);
850 CHKERR DMMoFEMTSSetIJacobian(dm, DM_NO_ELEMENT, null, schurAssembly, null);
851 CHKERR DMMoFEMTSSetIJacobian(dm, elementVolumeName, elasticFeLhs, null,
852 null);
853 CHKERR DMMoFEMTSSetIJacobian(dm, naturalBcElement, elasticBcLhs, null,
854 null);
855 CHKERR DMMoFEMTSSetIJacobian(dm, DM_NO_ELEMENT, null, null, schurAssembly);
856 }
857 MoFEMFunctionReturn(0);
858}
859
860MoFEMErrorCode EshelbianCore::setUpTSElastic(TS ts, Mat m, Vec f) {
861 MoFEMFunctionBegin;
862 boost::shared_ptr<TsCtx> ts_ctx;
863 CHKERR DMMoFEMGetTsCtx(dmElastic, ts_ctx);
864
865 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
866 CHKERR TSSetI2Function(ts, f, PETSC_NULL, PETSC_NULL);
867 CHKERR TSSetI2Jacobian(ts, m, m, PETSC_NULL, PETSC_NULL);
868 } else {
869 CHKERR TSSetIFunction(ts, f, PETSC_NULL, PETSC_NULL);
870 CHKERR TSSetIJacobian(ts, m, m, PETSC_NULL, PETSC_NULL);
871 }
872
873 CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx.get(), PETSC_NULL);
874
875 auto add_schur_streach_op = [this](auto &list, SmartPetscObj<Mat> S,
876 SmartPetscObj<AO> aoS) {
877 MoFEMFunctionBegin;
878 for (auto &fe : list)
879 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
880 fe_cast->addStreachSchurMatrix(S, aoS);
881 else
882 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
883 MoFEMFunctionReturn(0);
884 };
885
886 auto add_schur_streach_pre = [this](auto &list, SmartPetscObj<Mat> S,
887 SmartPetscObj<AO> aoS) {
888 MoFEMFunctionBegin;
889 for (auto &fe : list)
890 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
891 fe_cast->addStreachSchurMatrix(S, aoS);
892 else
893 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
894 MoFEMFunctionReturn(0);
895 };
896
897 auto add_schur_bubble_op = [this](auto &list, SmartPetscObj<Mat> S,
898 SmartPetscObj<AO> aoS) {
899 MoFEMFunctionBegin;
900 for (auto &fe : list)
901 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
902 fe_cast->addBubbleSchurMatrix(S, aoS);
903 else
904 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
905 MoFEMFunctionReturn(0);
906 };
907
908 auto add_schur_bubble_pre = [this](auto &list, SmartPetscObj<Mat> S,
909 SmartPetscObj<AO> aoS) {
910 MoFEMFunctionBegin;
911 for (auto &fe : list)
912 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
913 fe_cast->addBubbleSchurMatrix(S, aoS);
914 else
915 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
916 MoFEMFunctionReturn(0);
917 };
918
919 auto add_schur_spatial_disp_op = [this](auto &list, SmartPetscObj<Mat> S,
920 SmartPetscObj<AO> aoS) {
921 MoFEMFunctionBegin;
922 for (auto &fe : list)
923 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
924 fe_cast->addSpatialDispStressSchurMatrix(S, aoS);
925 else
926 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
927 MoFEMFunctionReturn(0);
928 };
929
930 auto add_schur_spatial_disp_pre = [this](auto &list, SmartPetscObj<Mat> S,
931 SmartPetscObj<AO> aoS) {
932 MoFEMFunctionBegin;
933 for (auto &fe : list)
934 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
935 fe_cast->addSpatialDispStressSchurMatrix(S, aoS);
936 else
937 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
938 MoFEMFunctionReturn(0);
939 };
940
941 auto add_schur_omega_op = [this](auto &list, SmartPetscObj<Mat> S,
942 SmartPetscObj<AO> aoS) {
943 MoFEMFunctionBegin;
944 for (auto &fe : list)
945 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
946 fe_cast->addOmegaSchurMatrix(S, aoS);
947 else
948 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
949 MoFEMFunctionReturn(0);
950 };
951
952 auto add_schur_omega_pre = [this](auto &list, SmartPetscObj<Mat> S,
953 SmartPetscObj<AO> ao) {
954 MoFEMFunctionBegin;
955 for (auto &fe : list)
956 if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
957 fe_cast->addOmegaSchurMatrix(S, ao);
958 else
959 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
960 MoFEMFunctionReturn(0);
961 };
962
963 const MoFEM::Problem *schur_streach_prb_ptr;
964 CHKERR DMMoFEMGetProblemPtr(dmElasticSchurStreach, &schur_streach_prb_ptr);
965 if (auto sub_data = schur_streach_prb_ptr->subProblemData) {
966 SmartPetscObj<Mat> Suu;
967 CHKERR DMCreateMatrix_MoFEM(dmElasticSchurStreach, Suu);
968 SmartPetscObj<AO> aoSuu = sub_data->getSmartRowMap();
969
970 CHKERR add_schur_streach_op(ts_ctx->loopsIJacobian, Suu, aoSuu);
971 CHKERR add_schur_streach_pre(ts_ctx->preProcessIJacobian, Suu, aoSuu);
972 CHKERR add_schur_streach_pre(ts_ctx->postProcessIJacobian, Suu, aoSuu);
973
974 const MoFEM::Problem *schur_bubble_prb_ptr;
975 CHKERR DMMoFEMGetProblemPtr(dmElasticSchurBubble, &schur_bubble_prb_ptr);
976 if (auto bubble_data = schur_bubble_prb_ptr->subProblemData) {
977 SmartPetscObj<Mat> SBubble;
978 CHKERR DMCreateMatrix_MoFEM(dmElasticSchurBubble, SBubble);
979 SmartPetscObj<AO> aoSBubble = bubble_data->getSmartRowMap();
980 CHKERR add_schur_bubble_op(ts_ctx->loopsIJacobian, SBubble,
981 aoSBubble);
982 CHKERR add_schur_bubble_pre(ts_ctx->preProcessIJacobian, SBubble,
983 aoSBubble);
984 CHKERR add_schur_bubble_pre(ts_ctx->postProcessIJacobian, SBubble,
985 aoSBubble);
986
987 const MoFEM::Problem *schur_omega_prb_ptr;
988 CHKERR DMMoFEMGetProblemPtr(dmElasticSchurOmega, &schur_omega_prb_ptr);
989 if (auto tet_stress_data = schur_omega_prb_ptr->subProblemData) {
990 SmartPetscObj<Mat> SOmega;
991 CHKERR DMCreateMatrix_MoFEM(dmElasticSchurOmega, SOmega);
992 SmartPetscObj<AO> aoSOmega = tet_stress_data->getSmartRowMap();
993 CHKERR add_schur_omega_op(ts_ctx->loopsIJacobian, SOmega,
994 aoSOmega);
995 CHKERR add_schur_omega_pre(ts_ctx->preProcessIJacobian, SOmega,
996 aoSOmega);
997 CHKERR add_schur_omega_pre(ts_ctx->postProcessIJacobian, SOmega,
998 aoSOmega);
999
1000 const MoFEM::Problem *schur_spatial_disp_prb_ptr;
1001 CHKERR DMMoFEMGetProblemPtr(dmElasticSchurSpatialDisp,
1002 &schur_spatial_disp_prb_ptr);
1003 if (auto spatial_disp_data =
1004 schur_spatial_disp_prb_ptr->subProblemData) {
1005 SmartPetscObj<Mat> Sw;
1006 CHKERR DMCreateMatrix_MoFEM(dmElasticSchurSpatialDisp, Sw);
1007 SmartPetscObj<AO> aoSw = spatial_disp_data->getSmartRowMap();
1008 CHKERR add_schur_spatial_disp_op(ts_ctx->loopsIJacobian, Sw,
1009 aoSw);
1010 CHKERR add_schur_spatial_disp_pre(ts_ctx->preProcessIJacobian, Sw,
1011 aoSw);
1012 CHKERR add_schur_spatial_disp_pre(ts_ctx->postProcessIJacobian, Sw,
1013 aoSw);
1014 } else
1015 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
1016 "Problem does not have sub-problem data");
1017
1018 } else
1019 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
1020 "Problem does not have sub-problem data");
1021
1022 } else
1023 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
1024 "Problem does not have sub-problem data");
1025
1026 } else
1027 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
1028 "Problem does not have sub-problem data");
1029
1030 struct Monitor : public FEMethod {
1031
1032 using Ele = ForcesAndSourcesCore;
1033 using VolEle = VolumeElementForcesAndSourcesCore;
1034 using VolOp = VolumeElementForcesAndSourcesCore::UserDataOperator;
1035 using SetPtsData = FieldEvaluatorInterface::SetPtsData;
1036
1037 EshelbianCore &eP;
1038 boost::shared_ptr<SetPtsData> dataFieldEval;
1039 boost::shared_ptr<VolEle> volPostProcEnergy;
1040 boost::shared_ptr<double> gEnergy;
1041
1042 Monitor(EshelbianCore &ep)
1043 : eP(ep),
1044 dataFieldEval(ep.mField.getInterface<FieldEvaluatorInterface>()
1045 ->getData<VolEle>()),
1046 volPostProcEnergy(new VolEle(ep.mField)), gEnergy(new double) {
1047 ierr = ep.mField.getInterface<FieldEvaluatorInterface>()->buildTree3D(
1048 dataFieldEval, "EP");
1049 CHKERRABORT(PETSC_COMM_WORLD, ierr);
1050
1051 auto no_rule = [](int, int, int) { return -1; };
1052
1053 auto set_element_for_field_eval = [&]() {
1054 boost::shared_ptr<Ele> vol_ele(dataFieldEval->feMethodPtr.lock());
1055 vol_ele->getRuleHook = no_rule;
1056 vol_ele->getUserPolynomialBase() =
1057 boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
1058 vol_ele->getOpPtrVector().push_back(new OpL2Transform());
1059
1060 vol_ele->getOpPtrVector().push_back(
1061 new OpCalculateHVecTensorField<3, 3>(
1062 eP.piolaStress, eP.dataAtPts->getApproxPAtPts()));
1063 vol_ele->getOpPtrVector().push_back(
1064 new OpCalculateHTensorTensorField<3, 3>(
1065 eP.bubbleField, eP.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
1066 vol_ele->getOpPtrVector().push_back(
1067 new OpCalculateTensor2SymmetricFieldValues<3>(
1068 eP.streachTensor, eP.dataAtPts->getLogStreachTensorAtPts(),
1069 MBTET));
1070 vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1071 eP.rotAxis, eP.dataAtPts->getRotAxisAtPts(), MBTET));
1072 vol_ele->getOpPtrVector().push_back(
1073 new OpCalculateTensor2FieldValues<3, 3>(
1074 eP.materialGradient, eP.dataAtPts->getBigGAtPts(), MBTET));
1075 vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1076 eP.spatialDisp, eP.dataAtPts->getSmallWAtPts(), MBTET));
1077 vol_ele->getOpPtrVector().push_back(
1078 new OpCalculateRotationAndSpatialGradient(eP.rotAxis,
1079 eP.dataAtPts));
1080 };
1081
1082 auto set_element_for_post_process = [&]() {
1083 volPostProcEnergy->getRuleHook = VolRule();
1084 volPostProcEnergy->getUserPolynomialBase() =
1085 boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
1086 volPostProcEnergy->getOpPtrVector().push_back(new OpL2Transform());
1087
1088 volPostProcEnergy->getOpPtrVector().push_back(
1089 new OpCalculateHVecTensorField<3, 3>(
1090 eP.piolaStress, eP.dataAtPts->getApproxPAtPts()));
1091 volPostProcEnergy->getOpPtrVector().push_back(
1092 new OpCalculateHTensorTensorField<3, 3>(
1093 eP.bubbleField, eP.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
1094 volPostProcEnergy->getOpPtrVector().push_back(
1095 new OpCalculateTensor2SymmetricFieldValues<3>(
1096 eP.streachTensor, eP.dataAtPts->getLogStreachTensorAtPts(),
1097 MBTET));
1098 volPostProcEnergy->getOpPtrVector().push_back(
1099 new OpCalculateVectorFieldValues<3>(
1100 eP.rotAxis, eP.dataAtPts->getRotAxisAtPts(), MBTET));
1101 volPostProcEnergy->getOpPtrVector().push_back(
1102 new OpCalculateTensor2FieldValues<3, 3>(
1103 eP.materialGradient, eP.dataAtPts->getBigGAtPts(), MBTET));
1104 volPostProcEnergy->getOpPtrVector().push_back(
1105 new OpCalculateVectorFieldValues<3>(
1106 eP.spatialDisp, eP.dataAtPts->getSmallWAtPts(), MBTET));
1107 volPostProcEnergy->getOpPtrVector().push_back(
1108 new OpCalculateRotationAndSpatialGradient(eP.rotAxis,
1109 eP.dataAtPts));
1110 volPostProcEnergy->getOpPtrVector().push_back(
1111 new OpCalculateStrainEnergy(eP.spatialDisp, eP.dataAtPts, gEnergy));
1112 };
1113
1114 set_element_for_field_eval();
1115 set_element_for_post_process();
1116 }
1117
1118 MoFEMErrorCode preProcess() { return 0; }
1119
1120 MoFEMErrorCode operator()() { return 0; }
1121
1122 MoFEMErrorCode postProcess() {
1123 MoFEMFunctionBegin;
1124
1125 // auto get_str_time = [](auto ts_t) {
1126 // std::ostringstream ss;
1127 // ss << boost::str(boost::format("%d") %
1128 // static_cast<int>(std::ceil(ts_t * 1e6)));
1129 // std::string s = ss.str();
1130 // return s;
1131 // };
1132
1133 auto get_step = [](auto ts_step) {
1134 std::ostringstream ss;
1135 ss << boost::str(boost::format("%d") % static_cast<int>(ts_step));
1136 std::string s = ss.str();
1137 return s;
1138 };
1139
1140 PetscViewer viewer;
1141 CHKERR PetscViewerBinaryOpen(
1142 PETSC_COMM_WORLD, ("restart_" + get_step(ts_step) + ".dat").c_str(),
1143 FILE_MODE_WRITE, &viewer);
1144 CHKERR VecView(ts_u, viewer);
1145 CHKERR PetscViewerDestroy(&viewer);
1146
1147 CHKERR eP.postProcessResults(1, "out_sol_elastic_" + get_step(ts_step) +
1148 ".h5m");
1149
1150 // Loop boundary elements with traction boundary conditions
1151 *gEnergy = 0;
1152 CHKERR eP.mField.loop_finite_elements(problemPtr->getName(), "EP",
1153 *volPostProcEnergy);
1154
1155 double body_energy;
1156 MPI_Allreduce(gEnergy.get(), &body_energy, 1, MPI_DOUBLE, MPI_SUM,
1157 eP.mField.get_comm());
1158 MOFEM_LOG_C("EP", Sev::inform, "Step %d time %3.4g strain energy %3.6e",
1159 ts_step, ts_t, body_energy);
1160
1161 auto post_proc_at_points = [&](std::array<double, 3> point,
1162 std::string str) {
1163 MoFEMFunctionBegin;
1164
1165 dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
1166
1167 struct OpPrint : public VolOp {
1168
1169 EshelbianCore &eP;
1170 std::array<double, 3> point;
1171 std::string str;
1172
1173 OpPrint(EshelbianCore &ep, std::array<double, 3> &point,
1174 std::string &str)
1175 : VolOp(ep.spatialDisp, VolOp::OPROW), eP(ep), point(point),
1176 str(str) {}
1177
1178 MoFEMErrorCode doWork(int side, EntityType type,
1179 DataForcesAndSourcesCore::EntData &data) {
1180 MoFEMFunctionBegin;
1181 if (type == MBTET) {
1182 if (getGaussPts().size2()) {
1183
1184 auto t_h = getFTensor2FromMat<3, 3>(eP.dataAtPts->hAtPts);
1185 auto t_approx_P =
1186 getFTensor2FromMat<3, 3>(eP.dataAtPts->approxPAtPts);
1187
1188 FTensor::Index<'i', 3> i;
1189 FTensor::Index<'j', 3> j;
1190 FTensor::Index<'k', 3> k;
1191 const double jac = determinantTensor3by3(t_h);
1192 FTensor::Tensor2<double, 3, 3> t_cauchy;
1193 t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
1194
1195 auto add = [&]() {
1196 std::ostringstream s;
1197 s << str << " elem " << getFEEntityHandle() << " ";
1198 return s.str();
1199 };
1200
1201 auto print_tensor = [](auto &t) {
1202 std::ostringstream s;
1203 s << t;
1204 return s.str();
1205 };
1206
1207 std::ostringstream print;
1208 MOFEM_LOG("EPSYNC", Sev::inform)
1209 << add() << "comm rank " << eP.mField.get_comm_rank();
1210 MOFEM_LOG("EPSYNC", Sev::inform)
1211 << add() << "point " << getVectorAdaptor(point.data(), 3);
1212 MOFEM_LOG("EPSYNC", Sev::inform)
1213 << add() << "coords at gauss pts " << getCoordsAtGaussPts();
1214 MOFEM_LOG("EPSYNC", Sev::inform)
1215 << add() << "w " << eP.dataAtPts->wAtPts;
1216 MOFEM_LOG("EPSYNC", Sev::inform)
1217 << add() << "Piola " << eP.dataAtPts->approxPAtPts;
1218 MOFEM_LOG("EPSYNC", Sev::inform)
1219 << add() << "Cauchy " << print_tensor(t_cauchy);
1220 }
1221 }
1222 MoFEMFunctionReturn(0);
1223 }
1224 };
1225
1226 if (auto fe_ptr = dataFieldEval->feMethodPtr.lock()) {
1227
1228 fe_ptr->getOpPtrVector().push_back(new OpPrint(eP, point, str));
1229 CHKERR eP.mField.getInterface<FieldEvaluatorInterface>()
1230 ->evalFEAtThePoint3D(
1231 point.data(), 1e-12, problemPtr->getName(), "EP",
1232 dataFieldEval, eP.mField.get_comm_rank(),
1233 eP.mField.get_comm_rank(), nullptr, MF_EXIST, QUIET);
1234 fe_ptr->getOpPtrVector().pop_back();
1235 }
1236
1237 MoFEMFunctionReturn(0);
1238 };
1239
1240 // Points for Cook beam
1241 std::array<double, 3> pointA = {48., 60., 5.};
1242 CHKERR post_proc_at_points(pointA, "Point A");
1243 MOFEM_LOG_SYNCHRONISE(eP.mField.get_comm());
1244
1245 std::array<double, 3> pointB = {48. / 2., 44. + (60. - 44.) / 2., 0.};
1246 CHKERR post_proc_at_points(pointB, "Point B");
1247 MOFEM_LOG_SYNCHRONISE(eP.mField.get_comm());
1248
1249 std::array<double, 3> pointC = {48. / 2., (44. - 0.) / 2., 0.};
1250 CHKERR post_proc_at_points(pointC, "Point C");
1251 MOFEM_LOG_SYNCHRONISE(eP.mField.get_comm());
1252
1253 MoFEMFunctionReturn(0);
1254 }
1255 };
1256
1257 boost::shared_ptr<FEMethod> monitor_ptr(new Monitor(*this));
1258 ts_ctx->getLoopsMonitor().push_back(
1259 TsCtx::PairNameFEMethodPtr(elementVolumeName, monitor_ptr));
1260
1261 CHKERR TSAppendOptionsPrefix(ts, "elastic_");
1262 CHKERR TSSetFromOptions(ts);
1263 MoFEMFunctionReturn(0);
1264}
1265
1266MoFEMErrorCode EshelbianCore::solveElastic(TS ts, Vec x) {
1267 MoFEMFunctionBegin;
1268
1269 CHKERR TSSetDM(ts, dmElastic);
1270
1271 SNES snes;
1272 CHKERR TSGetSNES(ts, &snes);
1273
1274 PetscViewerAndFormat *vf;
1275 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
1276 PETSC_VIEWER_DEFAULT, &vf);
1277 CHKERR SNESMonitorSet(
1278 snes,
1279 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
1280 vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
1281
1282 PetscSection section;
1283 CHKERR DMGetDefaultSection(dmElastic, &section);
1284 int num_fields;
1285 CHKERR PetscSectionGetNumFields(section, &num_fields);
1286 for (int ff = 0; ff != num_fields; ff++) {
1287 const char *field_name;
1288 CHKERR PetscSectionGetFieldName(section, ff, &field_name);
1289 MOFEM_LOG_C("EP", Sev::inform, "Field %d name %s", ff, field_name);
1290 }
1291
1292 CHKERR DMoFEMMeshToLocalVector(dmElastic, x, INSERT_VALUES, SCATTER_FORWARD);
1293
1294 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
1295 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
1296
1297
1298
1299 // Adding field split solver
1300
1301 KSP ksp;
1302 CHKERR SNESGetKSP(snes, &ksp);
1303 PC pc;
1304 CHKERR KSPGetPC(ksp, &pc);
1305 PetscBool is_uu_field_split;
1306 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_uu_field_split);
1307 if (is_uu_field_split) {
1308
1309 const MoFEM::Problem *schur_uu_ptr;
1310 CHKERR DMMoFEMGetProblemPtr(dmElasticSchurStreach, &schur_uu_ptr);
1311 if (auto uu_data = schur_uu_ptr->subProblemData) {
1312
1313 const MoFEM::Problem *prb_ptr;
1314 CHKERR DMMoFEMGetProblemPtr(dmElastic, &prb_ptr);
1315 map<std::string, IS> is_map;
1316 for (int ff = 0; ff != num_fields; ff++) {
1317 const char *field_name;
1318 CHKERR PetscSectionGetFieldName(section, ff, &field_name);
1319 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1320 prb_ptr->getName(), ROW, field_name, 0, MAX_DOFS_ON_ENTITY,
1321 &is_map[field_name]);
1322 }
1323 // CHKERR mField.getInterface<ISManager>()
1324 // ->isCreateProblemFieldAndEntityType(
1325 // prb_ptr->getName(), ROW, piolaStress, MBTET, MBTET, 0,
1326 // MAX_DOFS_ON_ENTITY, &is_map["T_STRESS_ON_TETS"]);
1327 // CHKERR mField.getInterface<ISManager>()
1328 // ->isCreateProblemFieldAndEntityType(
1329 // prb_ptr->getName(), ROW, piolaStress, MBTRI, MBTRI, 0,
1330 // MAX_DOFS_ON_ENTITY, &is_map["T_STRESS_ON_TRIS"]);
1331
1332 CHKERR uu_data->getRowIs(&is_map["E_IS_SUU"]);
1333
1334 CHKERR PCFieldSplitSetIS(pc, NULL, is_map[streachTensor]);
1335 CHKERR PCFieldSplitSetIS(pc, NULL, is_map["E_IS_SUU"]);
1336
1337 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER,
1338 schurAssembly->Suu);
1339
1340 CHKERR PCSetUp(pc);
1341 PetscInt n;
1342 KSP *uu_ksp;
1343 CHKERR PCFieldSplitGetSubKSP(pc, &n, &uu_ksp);
1344 PC bubble_pc;
1345 CHKERR KSPGetPC(uu_ksp[1], &bubble_pc);
1346 PetscBool is_bubble_field_split;
1347 PetscObjectTypeCompare((PetscObject)bubble_pc, PCFIELDSPLIT,
1348 &is_bubble_field_split);
1349 if (is_bubble_field_split) {
1350
1351 const MoFEM::Problem *schur_bubble_ptr;
1352 CHKERR DMMoFEMGetProblemPtr(dmElasticSchurBubble, &schur_bubble_ptr);
1353 if (auto bubble_data = schur_bubble_ptr->subProblemData) {
1354
1355 CHKERR bubble_data->getRowIs(&is_map["E_IS_BUBBLE"]);
1356
1357 AO uu_ao;
1358 CHKERR uu_data->getRowMap(&uu_ao);
1359
1360 CHKERR AOApplicationToPetscIS(uu_ao, is_map[bubbleField]);
1361 CHKERR PCFieldSplitSetIS(bubble_pc, NULL, is_map[bubbleField]);
1362 CHKERR PCFieldSplitSetIS(bubble_pc, NULL, is_map["E_IS_BUBBLE"]);
1363 CHKERR PCFieldSplitSetSchurPre(
1364 bubble_pc, PC_FIELDSPLIT_SCHUR_PRE_USER, schurAssembly->SBubble);
1365
1366 CHKERR PCSetUp(bubble_pc);
1367 PetscInt bubble_n;
1368 KSP *bubble_ksp;
1369 CHKERR PCFieldSplitGetSubKSP(bubble_pc, &bubble_n, &bubble_ksp);
1370 PC omega_pc;
1371 CHKERR KSPGetPC(bubble_ksp[1], &omega_pc);
1372 PetscBool is_omega_field_split;
1373 PetscObjectTypeCompare((PetscObject)omega_pc, PCFIELDSPLIT,
1374 &is_omega_field_split);
1375
1376 if (is_omega_field_split) {
1377
1378 const MoFEM::Problem *schur_omega_ptr;
1379 CHKERR DMMoFEMGetProblemPtr(dmElasticSchurOmega, &schur_omega_ptr);
1380 if (auto omega_data = schur_omega_ptr->subProblemData) {
1381
1382 AO bubble_ao;
1383 CHKERR bubble_data->getRowMap(&bubble_ao);
1384
1385 CHKERR AOApplicationToPetscIS(uu_ao, is_map[rotAxis]);
1386 CHKERR AOApplicationToPetscIS(bubble_ao, is_map[rotAxis]);
1387 CHKERR omega_data->getRowIs(&is_map["E_IS_OMEGA"]);
1388
1389 CHKERR PCFieldSplitSetIS(omega_pc, NULL, is_map[rotAxis]);
1390 CHKERR PCFieldSplitSetIS(omega_pc, NULL, is_map["E_IS_OMEGA"]);
1391 CHKERR PCFieldSplitSetSchurPre(omega_pc,
1392 PC_FIELDSPLIT_SCHUR_PRE_USER,
1393 schurAssembly->SOmega);
1394
1395 CHKERR PCSetUp(omega_pc);
1396 PetscInt omega_n;
1397 KSP *omega_ksp;
1398 CHKERR PCFieldSplitGetSubKSP(omega_pc, &omega_n, &omega_ksp);
1399 PC w_pc;
1400 CHKERR KSPGetPC(omega_ksp[1], &w_pc);
1401 PetscBool is_w_field_split;
1402 PetscObjectTypeCompare((PetscObject)w_pc, PCFIELDSPLIT,
1403 &is_w_field_split);
1404 if (is_w_field_split) {
1405
1406 const MoFEM::Problem *schur_w_ptr;
1407 CHKERR DMMoFEMGetProblemPtr(dmElasticSchurSpatialDisp,
1408 &schur_w_ptr);
1409 if (auto w_data = schur_w_ptr->subProblemData) {
1410
1411 AO omega_ao;
1412 CHKERR omega_data->getRowMap(&omega_ao);
1413
1414 CHKERR AOApplicationToPetscIS(uu_ao, is_map[spatialDisp]);
1415 CHKERR AOApplicationToPetscIS(bubble_ao, is_map[spatialDisp]);
1416 CHKERR AOApplicationToPetscIS(omega_ao, is_map[spatialDisp]);
1417 CHKERR w_data->getRowIs(&is_map["E_IS_W"]);
1418
1419 CHKERR PCFieldSplitSetIS(w_pc, NULL, is_map[spatialDisp]);
1420 CHKERR PCFieldSplitSetIS(w_pc, NULL, is_map["E_IS_W"]);
1421 CHKERR PCFieldSplitSetSchurPre(
1422 w_pc, PC_FIELDSPLIT_SCHUR_PRE_USER, schurAssembly->Sw);
1423
1424 CHKERR AODestroy(&omega_ao);
1425 }
1426 }
1427
1428 CHKERR PetscFree(omega_ksp);
1429 }
1430 }
1431 CHKERR PetscFree(bubble_ksp);
1432 CHKERR AODestroy(&uu_ao);
1433 }
1434 }
1435 CHKERR PetscFree(uu_ksp);
1436
1437 for (auto &m : is_map)
1438 CHKERR ISDestroy(&m.second);
1439 }
1440 }
1441
1442 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
1443 Vec xx;
1444 CHKERR VecDuplicate(x, &xx);
1445 CHKERR VecZeroEntries(xx);
1446 CHKERR TS2SetSolution(ts, x, xx);
1447 CHKERR VecDestroy(&xx);
1448 } else {
1449 CHKERR TSSetSolution(ts, x);
1450 }
1451
1452 CHKERR TSSolve(ts, PETSC_NULL);
1453
1454 // CHKERR TSGetSNES(ts, &snes);
1455 int lin_solver_iterations;
1456 CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
1457 MOFEM_LOG("EP", Sev::inform)
1458 << "Number of linear solver iterations " << lin_solver_iterations;
1459
1460 PetscBool test_cook_flg = PETSC_FALSE;
1461 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test_cook", &test_cook_flg,
1462 PETSC_NULL);
1463 if (test_cook_flg)
1464 if (lin_solver_iterations != 2)
1465 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1466 "Expected number of iterations is different than expected");
1467
1468 MoFEMFunctionReturn(0);
1469}
1470
1471MoFEMErrorCode EshelbianCore::postProcessResults(const int tag,
1472 const std::string file) {
1473 MoFEMFunctionBegin;
1474
1475 if (!dataAtPts) {
1476 dataAtPts =
1477 boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
1478 }
1479
1480 PostProcVolumeOnRefinedMesh post_proc(mField);
1481 post_proc.getUserPolynomialBase() =
1482 boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
1483 post_proc.getOpPtrVector().push_back(new OpL2Transform());
1484
1485 CHKERR post_proc.generateReferenceElementMesh();
1486 post_proc.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
1487 piolaStress, dataAtPts->getApproxPAtPts()));
1488 post_proc.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
1489 bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
1490 post_proc.getOpPtrVector().push_back(
1491 new OpCalculateTensor2SymmetricFieldValues<3>(
1492 streachTensor, dataAtPts->getLogStreachTensorAtPts(), MBTET));
1493 post_proc.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1494 rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
1495 post_proc.getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
1496 materialGradient, dataAtPts->getBigGAtPts(), MBTET));
1497 post_proc.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1498 spatialDisp, dataAtPts->getSmallWAtPts(), MBTET));
1499
1500 // evaluate derived quantities
1501 post_proc.getOpPtrVector().push_back(
1502 new OpCalculateRotationAndSpatialGradient(rotAxis, dataAtPts));
1503
1504 // evaluate integration points
1505 post_proc.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
1506 spatialDisp, tag, true, false, dataAtPts, physicalEquations));
1507 post_proc.getOpPtrVector().push_back(new OpPostProcDataStructure(
1508 post_proc.postProcMesh, post_proc.mapGaussPts, spatialDisp, dataAtPts));
1509
1510 CHKERR DMoFEMLoopFiniteElements(dM, elementVolumeName.c_str(), &post_proc);
1511 CHKERR post_proc.writeFile(file.c_str());
1512 MoFEMFunctionReturn(0);
1513}
1514
1515} // namespace EshelbianPlasticity
p
static Index< 'p', 3 > p
Definition: BasicFeTools.hpp:90
CGGTonsorialBubbleBase.hpp
Implementation of tonsorial bubble base div(v) = 0.
NBVOLUMETET_CCG_BUBBLE
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
Definition: CGGTonsorialBubbleBase.hpp:19
DM_NO_ELEMENT
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:10
EshelbianPlasticity.hpp
Eshelbian plasticity interface.
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
QUIET
@ QUIET
Definition: definitions.h:208
NOISY
@ NOISY
Definition: definitions.h:211
ROW
@ ROW
Definition: definitions.h:123
MF_ZERO
@ MF_ZERO
Definition: definitions.h:98
MF_EXIST
@ MF_EXIST
Definition: definitions.h:100
MAX_DOFS_ON_ENTITY
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
USER_BASE
@ USER_BASE
user implemented approximation base
Definition: definitions.h:68
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
H1
@ H1
continuous field
Definition: definitions.h:85
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
dim
const int dim
Definition: elec_phys_2D.cpp:12
ShapeDiffMBTET
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:319
ShapeMBTET
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition: fem_tools.c:306
n
FTensor::Index< 'n', SPACE_DIM > n
Definition: fluid_structure_eigenproblem.cpp:66
m
FTensor::Index< 'm', SPACE_DIM > m
Definition: fluid_structure_eigenproblem.cpp:65
order
int order
Definition: fluid_structure_eigenproblem.cpp:84
phi
static double phi
Definition: fluid_structure_eigenproblem.cpp:69
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1111
MoFEM::DMMoFEMCreateSubDM
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMoFEM.cpp:221
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:485
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:444
MoFEM::DMMoFEMTSSetIFunction
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
Definition: DMMoFEM.cpp:788
MoFEM::DMMoFEMCreateMoFEM
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMoFEM.cpp:118
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:511
MoFEM::DMMoFEMGetProblemPtr
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:414
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:244
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1185
MoFEM::DMMoFEMGetTsCtx
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMoFEM.cpp:1130
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:574
MoFEM::DMMoFEMTSSetIJacobian
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
Definition: DMMoFEM.cpp:841
MoFEM::DMMoFEMTSSetI2Jacobian
PetscErrorCode DMMoFEMTSSetI2Jacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
Definition: DMMoFEM.cpp:1005
MoFEM::DMMoFEMTSSetI2Function
PetscErrorCode DMMoFEMTSSetI2Function(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS implicit function evaluation function
Definition: DMMoFEM.cpp:963
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::CoreInterface::modify_finite_element_add_field_col
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
MoFEM::CoreInterface::modify_finite_element_add_field_row
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
MoFEM::CoreInterface::get_field_structure
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::CoreInterface::set_field_order
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
MoFEM::CoreInterface::add_ents_to_field_by_type
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
PostProcTemplateOnRefineMesh::writeFile
MoFEMErrorCode writeFile(const std::string file_name, const char *file_type="MOAB", const char *file_options="PARALLEL=WRITE_PART")
wrote results in (MOAB) format, use "file_name.h5m"
Definition: PostProcOnRefMesh.hpp:253
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::CoreInterface::loop_finite_elements
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Definition: MeshsetsManager.hpp:71
MoFEM::ProblemsManager::removeDofsOnEntities
MoFEMErrorCode removeDofsOnEntities(const std::string problem_name, const std::string field_name, const Range ents, const int lo_coeff=0, const int hi_coeff=MAX_DOFS_ON_ENTITY, const int lo_order=0, const int hi_order=100, int verb=VERBOSE, const bool debug=false)
Remove DOFs from problem.
Definition: ProblemsManager.cpp:2607
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
ts_ctx
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1876
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
EshelbianPlasticity::EntData
DataForcesAndSourcesCore::EntData EntData
Definition: EshelbianPlasticity.hpp:28
EshelbianPlasticity::CGG_BubbleBase_MBTET
MoFEMErrorCode CGG_BubbleBase_MBTET(const int p, const double *N, const double *diffN, FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > &phi, const int gdim)
Calculate CGGT tonsorial bubble base.
Definition: CGGTonsorialBubbleBase.cpp:20
EshelbianPlasticity::f
double f(const double v)
Definition: EshelbianOperators.cpp:22
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::getVectorAdaptor
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:31
MoFEM::TsMonitorSet
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:219
MoFEM::DMMoFEMSetDestroyProblem
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMoFEM.cpp:426
MoFEM::DataForcesAndSourcesCore
DEPRECATED typedef EntitiesFieldData DataForcesAndSourcesCore
Definition: EntitiesFieldData.hpp:1760
MoFEM::createSmartDM
auto createSmartDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:139
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1352
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
EshelbianPlasticity::BcDisp::BcDisp
BcDisp(std::string name, std::vector< double > &attr, Range &faces)
Definition: EshelbianPlasticity.cpp:423
EshelbianPlasticity::BcRot::BcRot
BcRot(std::string name, std::vector< double > &attr, Range &faces)
Definition: EshelbianPlasticity.cpp:433
EshelbianPlasticity::CGGUserPolynomialBase::getValueHdivForCGGBubble
MoFEMErrorCode getValueHdivForCGGBubble(MatrixDouble &pts)
Definition: EshelbianPlasticity.cpp:587
EshelbianPlasticity::CGGUserPolynomialBase::getValue
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: EshelbianPlasticity.cpp:554
EshelbianPlasticity::CGGUserPolynomialBase::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, BaseFunctionUnknownInterface **iface) const
Definition: EshelbianPlasticity.cpp:548
EshelbianPlasticity::EshelbianCore::setUpTSElastic
MoFEMErrorCode setUpTSElastic(TS ts, Mat m, Vec f)
Definition: EshelbianPlasticity.cpp:860
EshelbianPlasticity::EshelbianCore::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Getting interface of core database.
Definition: EshelbianPlasticity.cpp:20
EshelbianPlasticity::EshelbianCore::dmElasticSchurBubble
SmartPetscObj< DM > dmElasticSchurBubble
Sub problem of dmElastic Schur.
Definition: EshelbianPlasticity.hpp:1182
EshelbianPlasticity::EshelbianCore::postProcessResults
MoFEMErrorCode postProcessResults(const int tag, const std::string file)
Definition: EshelbianPlasticity.cpp:1471
EshelbianPlasticity::EshelbianCore::addBoundaryFiniteElement
MoFEMErrorCode addBoundaryFiniteElement(const EntityHandle meshset=0)
Definition: EshelbianPlasticity.cpp:254
EshelbianPlasticity::EshelbianCore::elasticFeLhs
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > elasticFeLhs
Definition: EshelbianPlasticity.hpp:1174
EshelbianPlasticity::EshelbianCore::dmElasticSchurStreach
SmartPetscObj< DM > dmElasticSchurStreach
Sub problem of dmElastic Schur.
Definition: EshelbianPlasticity.hpp:1181
EshelbianPlasticity::EshelbianCore::addDMs
MoFEMErrorCode addDMs(const BitRefLevel bit=BitRefLevel().set(0))
Definition: EshelbianPlasticity.cpp:305
EshelbianPlasticity::EshelbianCore::setGenericFaceElementOps
MoFEMErrorCode setGenericFaceElementOps(const bool add_elastic, const bool add_material, boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > &fe_rhs, boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > &fe_lhs)
Definition: EshelbianPlasticity.cpp:781
EshelbianPlasticity::EshelbianCore::elasticBcRhs
boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > elasticBcRhs
Definition: EshelbianPlasticity.hpp:1176
EshelbianPlasticity::EshelbianCore::setElasticElementOps
MoFEMErrorCode setElasticElementOps(const int tag)
Definition: EshelbianPlasticity.cpp:811
EshelbianPlasticity::EshelbianCore::getTractionFreeBc
MoFEMErrorCode getTractionFreeBc(const EntityHandle meshset, boost::shared_ptr< TractionFreeBc > &bc_ptr, const std::string disp_block_set_name, const std::string rot_block_set_name)
Definition: EshelbianPlasticity.cpp:454
EshelbianPlasticity::EshelbianCore::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: EshelbianPlasticity.hpp:1170
EshelbianPlasticity::EshelbianCore::setGenericVolumeElementOps
MoFEMErrorCode setGenericVolumeElementOps(const int tag, const bool add_elastic, const bool add_material, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > &fe_rhs, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > &fe_lhs)
Definition: EshelbianPlasticity.cpp:697
EshelbianPlasticity::EshelbianCore::dmElasticSchurOmega
SmartPetscObj< DM > dmElasticSchurOmega
Sub problem of dmElastic Schur.
Definition: EshelbianPlasticity.hpp:1185
EshelbianPlasticity::EshelbianCore::bcSpatialRotationVecPtr
boost::shared_ptr< BcRotVec > bcSpatialRotationVecPtr
Definition: EshelbianPlasticity.hpp:1216
EshelbianPlasticity::EshelbianCore::dmElasticSchurSpatialDisp
SmartPetscObj< DM > dmElasticSchurSpatialDisp
Sub problem of dmElastic Schur.
Definition: EshelbianPlasticity.hpp:1184
EshelbianPlasticity::EshelbianCore::EshelbianCore
EshelbianCore(MoFEM::Interface &m_field)
Definition: EshelbianPlasticity.cpp:40
EshelbianPlasticity::EshelbianCore::dM
SmartPetscObj< DM > dM
Coupled problem all fields.
Definition: EshelbianPlasticity.hpp:1179
EshelbianPlasticity::EshelbianCore::solveElastic
MoFEMErrorCode solveElastic(TS ts, Vec x)
Definition: EshelbianPlasticity.cpp:1266
EshelbianPlasticity::EshelbianCore::addVolumeFiniteElement
MoFEMErrorCode addVolumeFiniteElement(const EntityHandle meshset=0)
Definition: EshelbianPlasticity.cpp:216
EshelbianPlasticity::EshelbianCore::setBaseVolumeElementOps
MoFEMErrorCode setBaseVolumeElementOps(const int tag, const bool do_rhs, const bool do_lhs, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > &fe)
Definition: EshelbianPlasticity.cpp:630
EshelbianPlasticity::EshelbianCore::dmElastic
SmartPetscObj< DM > dmElastic
Elastic problem.
Definition: EshelbianPlasticity.hpp:1180
EshelbianPlasticity::EshelbianCore::schurAssembly
boost::shared_ptr< EpFEMethod > schurAssembly
Definition: EshelbianPlasticity.hpp:1177
EshelbianPlasticity::EshelbianCore::bcSpatialDispVecPtr
boost::shared_ptr< BcDispVec > bcSpatialDispVecPtr
Definition: EshelbianPlasticity.hpp:1215
EshelbianPlasticity::EshelbianCore::physicalEquations
boost::shared_ptr< PhysicalEquations > physicalEquations
Definition: EshelbianPlasticity.hpp:1171
EshelbianPlasticity::EshelbianCore::elasticFeRhs
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > elasticFeRhs
Definition: EshelbianPlasticity.hpp:1173
EshelbianPlasticity::EshelbianCore::bcSpatialFreeTraction
boost::shared_ptr< TractionFreeBc > bcSpatialFreeTraction
Definition: EshelbianPlasticity.hpp:1218
EshelbianPlasticity::EshelbianCore::bcSpatialTraction
boost::shared_ptr< TractionBcVec > bcSpatialTraction
Definition: EshelbianPlasticity.hpp:1217
EshelbianPlasticity::EshelbianCore::addFields
MoFEMErrorCode addFields(const EntityHandle meshset=0)
Definition: EshelbianPlasticity.cpp:86
EshelbianPlasticity::EshelbianCore::elasticBcLhs
boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > elasticBcLhs
Definition: EshelbianPlasticity.hpp:1175
EshelbianPlasticity::FaceRule::operator()
int operator()(int p_row, int p_col, int p_data) const
Definition: EshelbianPlasticity.cpp:538
EshelbianPlasticity::OpJacobian::evaluateRhs
virtual MoFEMErrorCode evaluateRhs(EntData &data)=0
EshelbianPlasticity::OpJacobian::evaluateLhs
virtual MoFEMErrorCode evaluateLhs(EntData &data)=0
EshelbianPlasticity::OpJacobian::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: EshelbianPlasticity.cpp:26
EshelbianPlasticity::OpPostProcDataStructure
Definition: EshelbianPlasticity.hpp:1041
EshelbianPlasticity::OpSpatialConsistency_dBubble_domega
Definition: EshelbianPlasticity.hpp:1031
EshelbianPlasticity::OpSpatialConsistency_dP_domega
Definition: EshelbianPlasticity.hpp:1021
EshelbianPlasticity::OpSpatialEquilibrium_dw_dP
Definition: EshelbianPlasticity.hpp:896
EshelbianPlasticity::OpSpatialEquilibrium_dw_dw
Definition: EshelbianPlasticity.hpp:908
EshelbianPlasticity::OpSpatialPhysical_du_dBubble
Definition: EshelbianPlasticity.hpp:948
EshelbianPlasticity::OpSpatialPhysical_du_dP
Definition: EshelbianPlasticity.hpp:935
EshelbianPlasticity::OpSpatialPhysical_du_domega
Definition: EshelbianPlasticity.hpp:961
EshelbianPlasticity::OpSpatialPhysical_du_du
Definition: EshelbianPlasticity.hpp:922
EshelbianPlasticity::OpSpatialRotation_domega_dBubble
Definition: EshelbianPlasticity.hpp:997
EshelbianPlasticity::OpSpatialRotation_domega_dP
Definition: EshelbianPlasticity.hpp:985
EshelbianPlasticity::OpSpatialRotation_domega_domega
Definition: EshelbianPlasticity.hpp:1009
EshelbianPlasticity::TractionBc::TractionBc
TractionBc(std::string name, std::vector< double > &attr, Range &faces)
Definition: EshelbianPlasticity.cpp:442
EshelbianPlasticity::VolRule
Set integration rule on element.
Definition: EshelbianPlasticity.cpp:531
EshelbianPlasticity::VolRule::operator()
int operator()(int p_row, int p_col, int p_data) const
Definition: EshelbianPlasticity.cpp:532
MoFEM::BaseFunction
Base class if inherited used to calculate base functions.
Definition: BaseFunction.hpp:40
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::CoreInterface::check_finite_element
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreInterface::add_field
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::EntPolynomialBaseCtx
Class used to pass element data to calculate base functions on tet,triangle,edge.
Definition: EntPolynomialBaseCtx.hpp:22
MoFEM::EntPolynomialBaseCtx::bAse
const FieldApproximationBase bAse
Definition: EntPolynomialBaseCtx.hpp:38
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
MoFEM::FieldEvaluatorInterface
Field evaluator interface.
Definition: FieldEvaluator.hpp:21
MoFEM::FieldEvaluatorInterface::getData
boost::shared_ptr< SPD > getData(const double *ptr=nullptr, const int nb_eval_points=0, const double eps=1e-12, VERBOSITY_LEVELS verb=QUIET)
Get the Data object.
Definition: FieldEvaluator.hpp:107
MoFEM::Field
Provide data structure for (tensor) field approximation.
Definition: FieldMultiIndices.hpp:51
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
operator doWork function is executed on FE rows
Definition: ForcesAndSourcesCore.hpp:564
MoFEM::ForcesAndSourcesCore
structure to get information form mofem into EntitiesFieldData
Definition: ForcesAndSourcesCore.hpp:22
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::OpCalculateHTensorTensorField
Calculate tenor field using tensor base, i.e. Hdiv/Hcurl.
Definition: UserDataOperators.hpp:2574
MoFEM::OpCalculateHVecTensorDivergence
Calculate divergence of tonsorial field using vectorial base.
Definition: UserDataOperators.hpp:2634
MoFEM::OpCalculateHVecTensorField
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Definition: UserDataOperators.hpp:2515
MoFEM::OpCalculateTensor2FieldValues
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
Definition: UserDataOperators.hpp:871
MoFEM::OpCalculateTensor2SymmetricFieldValuesDot
Calculate symmetric tensor field rates ant integratio pts.
Definition: UserDataOperators.hpp:1068
MoFEM::OpCalculateTensor2SymmetricFieldValues
Calculate symmetric tensor field values at integration pts.
Definition: UserDataOperators.hpp:974
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:578
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:450
MoFEM::OpHOSetContravariantPiolaTransformOnFace3D
transform Hdiv base fluxes from reference element to physical triangle
Definition: HODataOperators.hpp:273
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:40
MoFEM::Problem::subProblemData
boost::shared_ptr< SubProblemData > subProblemData
Definition: ProblemsMultiIndices.hpp:91
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
MoFEM::SmartPetscObj
intrusive_ptr for managing petsc objects
Definition: PetscSmartObj.hpp:81
MoFEM::TsCtx::loopsIJacobian
FEMethodsSequence loopsIJacobian
Definition: TsCtx.hpp:27
MoFEM::TsCtx::getLoopsMonitor
FEMethodsSequence & getLoopsMonitor()
Get the loops to do Monitor object.
Definition: TsCtx.hpp:108
MoFEM::TsCtx::preProcessIJacobian
BasicMethodsSequence preProcessIJacobian
Definition: TsCtx.hpp:32
MoFEM::TsCtx::postProcessIJacobian
BasicMethodsSequence postProcessIJacobian
Definition: TsCtx.hpp:33
MoFEM::UnknownInterface
base class for all interface classes
Definition: UnknownInterface.hpp:34
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
Monitor
Monitor solution.
Definition: heat_equation.cpp:66
PostProcTemplateOnRefineMesh::mapGaussPts
std::vector< EntityHandle > mapGaussPts
Definition: PostProcOnRefMesh.hpp:125
PostProcTemplateOnRefineMesh::postProcMesh
moab::Interface & postProcMesh
Definition: PostProcOnRefMesh.hpp:122
PostProcTemplateVolumeOnRefinedMesh::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Definition: PostProcOnRefMesh.hpp:301
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:957
VolOp
VolEle::UserDataOperator VolOp
Definition: test_cache_on_entities.cpp:23
VolEle
VolumeElementForcesAndSourcesCore VolEle
Definition: test_cache_on_entities.cpp:22
Generated by Doxygen 1.9.5 and hosted at University of Glasgow