v0.16.3
Loading...
Searching...
No Matches
ElasticTie.hpp
Go to the documentation of this file.
1/**
2 * @file ElasticTie.hpp
3 *
4 * @brief TIE constraint support for the elastic tutorial
5 *
6 * @copyright Copyright (c) 2025
7 *
8 */
9
10#ifndef __ELASTIC_TIE_HPP__
11#define __ELASTIC_TIE_HPP__
12
13#include <ElasticExample.hpp>
14#include <regex>
15#include <quad.h>
16
17namespace ElasticTie {
18
19FTENSOR_INDEX(SPACE_DIM, i);
20
21// we know this because we created the prisms in insert_constraint_prisms.cpp
22constexpr int SLAVE_FACE_SIDE = 3;
23constexpr int MASTER_FACE_SIDE = 4;
24
25struct CommonData {
26 bool hasTieConstraints = false;
27 std::string tiePrismFeName = "TIE_PRISM_FE";
28 std::string dispFieldName = "U";
29 FieldApproximationBase dispFieldBase = AINSWORTH_LEGENDRE_BASE;
30 int dispFieldOrder = 1;
31 int tieRefLevel = 1;
32 boost::shared_ptr<Range> tieMasterFaces = nullptr;
33 boost::shared_ptr<Range> tieSlaveFaces = nullptr;
34 boost::shared_ptr<Range> tiePrisms = nullptr;
35 boost::shared_ptr<MatrixDouble> dispAtGaussPts =
36 boost::make_shared<MatrixDouble>();
37 boost::shared_ptr<MatrixDouble> lambdaAtGaussPts =
38 boost::make_shared<MatrixDouble>();
39 boost::shared_ptr<FlatPrismElementForcesAndSourcesCore> tiePrismRhsFE =
40 nullptr;
41 boost::shared_ptr<FlatPrismElementForcesAndSourcesCore> tiePrismLhsFE =
42 nullptr;
43};
44
45inline MoFEMErrorCode makeClosureMeshset(moab::Interface &moab,
46 const Range &ents, const int max_dim,
47 EntityHandle &meshset) {
48 MoFEMFunctionBegin;
49 Range closure(ents);
50 for (int d = 0; d < max_dim; ++d) {
51 Range adj;
52 CHKERR moab.get_adjacencies(ents, d, false, adj,
53 moab::Interface::UNION);
54 closure.merge(adj);
55 }
56 CHKERR moab.create_meshset(MESHSET_SET, meshset);
57 CHKERR moab.add_entities(meshset, closure);
58 MoFEMFunctionReturn(0);
59}
60
61inline Range getPhysicalVolumeEntities(const Range &ents) {
62 Range physical_ents = ents.subset_by_type(MBTET);
63 physical_ents.merge(ents.subset_by_type(MBHEX));
64 return physical_ents;
65}
66
67inline MoFEMErrorCode getPrismTriangleFaces(moab::Interface &moab,
68 const Range &ents,
69 Range &triangle_faces) {
70 MoFEMFunctionBegin;
71 triangle_faces.clear();
72 const auto prisms = ents.subset_by_type(MBPRISM);
73 if (!prisms.empty()) {
74 CHKERR moab.get_adjacencies(prisms, 2, false, triangle_faces,
75 moab::Interface::UNION);
76 triangle_faces = triangle_faces.subset_by_type(MBTRI);
77 }
78 MoFEMFunctionReturn(0);
79}
80
81/**
82 * @brief Load a mesh while preserving the shared entities needed by TIE
83 * constraints.
84 *
85 * @param mField MoFEM interface
86 * @param physicalBoundaryFaces Physical boundary faces owned by this process
87 */
88template <int DIM>
89inline MoFEMErrorCode loadTieMesh(MoFEM::Interface &mField,
90 Range &physicalBoundaryFaces) {
91 MoFEMFunctionBegin;
92 auto simple = mField.getInterface<Simple>();
93 auto meshset_mng = mField.getInterface<MeshsetsManager>();
94 auto comm_interface = mField.getInterface<CommInterface>();
95 CHKERR simple->getOptions();
96 physicalBoundaryFaces.clear();
97
98 if (DIM == 3 && mField.get_comm_size() > 1) {
99 auto load_fun = [&](MoFEM::Interface &mf, const char *file_name,
100 const char *) -> MoFEMErrorCode {
101 MoFEMFunctionBegin;
102 auto &moab = mf.get_moab();
103 CHKERR moab.load_file(file_name, 0,
104 "PARALLEL=BCAST;PARTITION=PARALLEL_PARTITION;");
105
106 auto *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
107 if (pcomm == nullptr)
108 pcomm = new ParallelComm(&moab, mf.get_comm());
109
110 const auto rank = pcomm->rank();
111 Tag part_tag = pcomm->part_tag();
112 Range all_ents, tagged_sets, proc_ents, off_proc_ents;
113 CHKERR moab.get_entities_by_handle(0, all_ents, false);
114 CHKERR moab.get_entities_by_type_and_tag(
115 0, MBENTITYSET, &part_tag, nullptr, 1, tagged_sets,
116 moab::Interface::UNION);
117
118 for (auto meshset : tagged_sets) {
119 int part = -1;
120 CHKERR moab.tag_get_data(part_tag, &meshset, 1, &part);
121 Range meshset_ents;
122 CHKERR moab.get_entities_by_handle(meshset, meshset_ents, true);
123 if (part == rank)
124 proc_ents.merge(meshset_ents);
125 else
126 off_proc_ents.merge(meshset_ents);
127 CHKERR moab.tag_clear_data(part_tag, meshset_ents, &part);
128 }
129
130 const auto all_volume_ents = all_ents.subset_by_dimension(3);
131 const auto proc_volume_ents = proc_ents.subset_by_dimension(3);
132 const auto all_phys_ents = getPhysicalVolumeEntities(all_ents);
133 const auto proc_phys_ents = getPhysicalVolumeEntities(proc_ents);
134
135 std::array<Range, 4> proc_ents_skin;
136 Skinner skin(&moab);
137 Range all_skin, proc_skin;
138 CHKERR skin.find_skin(0, all_volume_ents, false, all_skin);
139 CHKERR skin.find_skin(0, proc_volume_ents, false, proc_skin);
140 proc_ents_skin[2] = subtract(proc_skin, all_skin);
141
142 Range all_phys_skin, local_phys_faces, physical_boundary_faces;
143 CHKERR skin.find_skin(0, all_phys_ents.subset_by_dimension(3), false,
144 all_phys_skin);
145 CHKERR moab.get_adjacencies(proc_phys_ents, 2, false, local_phys_faces,
146 moab::Interface::UNION);
147 physical_boundary_faces = intersect(all_phys_skin, local_phys_faces);
148 if (!physical_boundary_faces.empty()) {
149 CHKERR moab.tag_clear_data(pcomm->partition_tag(),
150 physical_boundary_faces, &rank);
151 }
152 proc_ents_skin[2].merge(physical_boundary_faces);
153
154 Range prism_triangle_faces;
155 CHKERR getPrismTriangleFaces(moab, proc_ents, prism_triangle_faces);
156 proc_ents_skin[2].merge(prism_triangle_faces);
157 CHKERR moab.get_adjacencies(proc_ents_skin[2], 1, false,
158 proc_ents_skin[1], moab::Interface::UNION);
159 CHKERR moab.get_connectivity(proc_ents_skin[2], proc_ents_skin[0],
160 false);
161
162 auto to_remove = off_proc_ents;
163 for (int d = 2; d >= 0; --d)
164 to_remove = subtract(to_remove, proc_ents_skin[d]);
165
166 Range all_meshsets;
167 CHKERR moab.get_entities_by_type(0, MBENTITYSET, all_meshsets, true);
168 for (auto meshset : all_meshsets)
169 CHKERR moab.remove_entities(meshset, to_remove);
170
171 for (int d = 3; d > 0; --d) {
172 Range ents = to_remove.subset_by_dimension(d);
173 if (!ents.empty())
174 CHKERR moab.delete_entities(ents);
175 }
176
177 physicalBoundaryFaces = physical_boundary_faces;
178 CHKERR pcomm->resolve_shared_ents(0, proc_ents, 3, -1,
179 proc_ents_skin.data());
180
181 MoFEMFunctionReturn(0);
182 };
183
184 CHKERR simple->loadFile("", "", load_fun);
185 } else {
186 CHKERR simple->loadFile();
187 }
188
189 if (DIM == 3 && mField.get_comm_size() > 1)
190 CHKERR comm_interface->synchroniseEntities(physicalBoundaryFaces);
191
192 CHKERR meshset_mng->setMeshsetFromFile();
193 MoFEMFunctionReturn(0);
194}
195
196
197inline bool isInsideReferenceTriangle(const double xi, const double eta,
198 const double tol = 1e-8) {
199 return xi >= -tol && eta >= -tol && xi + eta <= 1 + tol;
200}
201
202inline bool isOnPrismTriFace(const int side, const EntityType type,
203 const int prism_face_side) {
204 switch (type) {
205 case MBVERTEX:
206 return prism_face_side == SLAVE_FACE_SIDE ? (side >= 0 && side < 3)
207 : (side >= 3 && side < 6);
208 case MBEDGE:
209 return prism_face_side == SLAVE_FACE_SIDE ? (side >= 0 && side < 3)
210 : (side >= 6 && side < 9);
211 case MBTRI:
212 return side == prism_face_side;
213 default:
214 return false;
215 }
216}
217
218inline bool isOnPrismTieTrace(const int side, const EntityType type) {
219 return isOnPrismTriFace(side, type, SLAVE_FACE_SIDE) ||
220 isOnPrismTriFace(side, type, MASTER_FACE_SIDE);
221}
222
223inline bool isOnSlaveTieRowBase(const int rr, const int row_side,
224 const EntityType row_type) {
225 switch (row_type) {
226 case MBVERTEX:
227 return rr < 3;
228 case MBEDGE:
229 return row_side >= 0 && row_side < 3;
230 case MBTRI:
231 return row_side == SLAVE_FACE_SIDE;
232 default:
233 return false;
234 }
235}
236
237struct SetIntegrationOnTiePrismFaces {
238
239 using FunRule = boost::function<int(int)>;
240 FunRule funRule = [](int p) { return 2 * p + 1; };
241 int tieRefLevel = 1;
242
243 SetIntegrationOnTiePrismFaces() = delete;
244
245 SetIntegrationOnTiePrismFaces(int tie_ref_level = 1)
246 : tieRefLevel(tie_ref_level) {}
247
248 SetIntegrationOnTiePrismFaces(FunRule fun_rule, int tie_ref_level = 1)
249 : funRule(fun_rule), tieRefLevel(tie_ref_level) {}
250
251 MoFEMErrorCode operator()(ForcesAndSourcesCore *fe_raw_ptr, int order_row,
252 int order_col, int order_data) {
253 MoFEMFunctionBegin;
254
255 auto fe_ptr = static_cast<FlatPrismElementForcesAndSourcesCore *>(fe_raw_ptr);
256 const int rule = funRule(std::max(order_data, std::max(order_row, order_col)));
257
258 if (tieRefLevel < 0) {
259 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
260 "-tie_ref_level has to be >= 0");
261 }
262
263 if (rule < 0 || rule >= QUAD_2D_TABLE_SIZE) {
264 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
265 "Triangle quadrature rule %d is not available for TIE prism face "
266 "integration",
267 rule);
268 }
269
270 if (QUAD_2D_TABLE[rule]->dim != 2) {
271 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
272 "Expected 2D quadrature for TIE prism face integration");
273 }
274
275 if (QUAD_2D_TABLE[rule]->order < rule) {
276 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
277 "Wrong quadrature order %d < %d for TIE prism face integration",
278 QUAD_2D_TABLE[rule]->order, rule);
279 }
280
281 CHKERR setBaseQuadrature(*fe_ptr, rule);
282
283 if (tieRefLevel > 0) {
284 const auto cache_key = std::make_pair(rule, tieRefLevel);
285 auto it = mapRefCoords.find(cache_key);
286 if (it == mapRefCoords.end()) {
287 MatrixDouble ref_gauss_pts;
288 CHKERR refineQuadrature(*fe_ptr, tieRefLevel, ref_gauss_pts);
289 mapRefCoords[cache_key].swap(ref_gauss_pts);
290 it = mapRefCoords.find(cache_key);
291 }
292 fe_ptr->gaussPts = it->second;
293 }
294
295 const MatrixDouble slave_gauss_pts = fe_ptr->gaussPts;
296 const auto nb_gauss_pts = slave_gauss_pts.size2();
297
298 auto &moab = fe_ptr->mField.get_moab();
299 auto get_prism_face_coords = [&](const int side, std::array<double, 9> &coords) {
300 MoFEMFunctionBegin;
301 const EntityHandle prism = fe_ptr->getFEEntityHandle();
302 const EntityHandle *prism_conn = nullptr;
303 int num_prism_nodes = 0;
304 CHKERR moab.get_connectivity(prism, prism_conn, num_prism_nodes, true);
305 if (num_prism_nodes != 6) {
306 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
307 "TIE prism is expected to have 6 nodes");
308 }
309
310 std::array<EntityHandle, 3> face_conn;
311 switch (side) {
312 case SLAVE_FACE_SIDE:
313 face_conn = {prism_conn[0], prism_conn[1], prism_conn[2]};
314 break;
315 case MASTER_FACE_SIDE:
316 face_conn = {prism_conn[3], prism_conn[4], prism_conn[5]};
317 break;
318 default:
319 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
320 "Unsupported TIE prism triangle side %d", side);
321 }
322
323 CHKERR moab.get_coords(face_conn.data(), 3, coords.data());
324 MoFEMFunctionReturn(0);
325 };
326
327 std::array<double, 9> slave_coords; // nodal coords
328 std::array<double, 9> master_coords; // nodal coords
329 CHKERR get_prism_face_coords(SLAVE_FACE_SIDE, slave_coords);
330 CHKERR get_prism_face_coords(MASTER_FACE_SIDE, master_coords);
331
332 MatrixDouble slave_global_coords(nb_gauss_pts, 3, false); // gp coords
333 MatrixDouble master_global_coords(nb_gauss_pts, 3, false); //gp coords
334 MatrixDouble slave_shape(nb_gauss_pts, 3, false);
335 CHKERR ShapeMBTRI(&slave_shape(0, 0), &slave_gauss_pts(0, 0),
336 &slave_gauss_pts(1, 0), nb_gauss_pts);
337 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
338 for (int dd = 0; dd != 3; ++dd) {
339 slave_global_coords(gg, dd) =
340 slave_shape(gg, 0) * slave_coords[0 + dd] +
341 slave_shape(gg, 1) * slave_coords[3 + dd] +
342 slave_shape(gg, 2) * slave_coords[6 + dd];
343 }
344 }
345
346 CHKERR projectSlavePointsToMasterPlane(master_coords.data(),
347 slave_global_coords,
348 master_global_coords);
349
350 MatrixDouble master_local_coords(nb_gauss_pts, 2, false);
351 CHKERR Tools::getLocalCoordinatesOnReferenceThreeNodeTri(
352 master_coords.data(), &master_global_coords(0, 0),
353 master_global_coords.size1(), &master_local_coords(0, 0));
354
355 MatrixDouble combined_gauss_pts(3, 2 * nb_gauss_pts, false);
356 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
357 const bool valid_master_point =
358 isInsideReferenceTriangle(master_local_coords(gg, 0),
359 master_local_coords(gg, 1));
360 const double paired_weight = valid_master_point ? slave_gauss_pts(2, gg) : 0.;
361
362 combined_gauss_pts(0, gg) = slave_gauss_pts(0, gg);
363 combined_gauss_pts(1, gg) = slave_gauss_pts(1, gg);
364 combined_gauss_pts(2, gg) = paired_weight;
365
366 combined_gauss_pts(0, gg + nb_gauss_pts) = master_local_coords(gg, 0);
367 combined_gauss_pts(1, gg + nb_gauss_pts) = master_local_coords(gg, 1);
368 combined_gauss_pts(2, gg + nb_gauss_pts) = paired_weight;
369 }
370 fe_ptr->gaussPts.swap(combined_gauss_pts);
371
372 MoFEMFunctionReturn(0);
373 }
374
375private:
376 static MoFEMErrorCode setBaseQuadrature(
377 FlatPrismElementForcesAndSourcesCore &fe, const int rule) {
378 MoFEMFunctionBegin;
379 const size_t nb_gauss_pts = QUAD_2D_TABLE[rule]->npoints;
380 fe.gaussPts.resize(3, nb_gauss_pts, false);
381 cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[1], 3,
382 &fe.gaussPts(0, 0), 1);
383 cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[2], 3,
384 &fe.gaussPts(1, 0), 1);
385 cblas_dcopy(nb_gauss_pts, QUAD_2D_TABLE[rule]->weights, 1,
386 &fe.gaussPts(2, 0), 1);
387 MoFEMFunctionReturn(0);
388 }
389
390 static MoFEMErrorCode projectPointToMasterPlane(const double *master_coords,
391 const double *point,
392 double *projected_point) {
393 MoFEMFunctionBegin;
394
395 // For the infinite master plane, the closest point is the orthogonal
396 // projection along the plane normal.
397 const double ax = master_coords[0];
398 const double ay = master_coords[1];
399 const double az = master_coords[2];
400 const double bx = master_coords[3];
401 const double by = master_coords[4];
402 const double bz = master_coords[5];
403 const double cx = master_coords[6];
404 const double cy = master_coords[7];
405 const double cz = master_coords[8];
406
407 const double abx = bx - ax;
408 const double aby = by - ay;
409 const double abz = bz - az;
410 const double acx = cx - ax;
411 const double acy = cy - ay;
412 const double acz = cz - az;
413
414 const double nx = aby * acz - abz * acy;
415 const double ny = abz * acx - abx * acz;
416 const double nz = abx * acy - aby * acx;
417 const double n_norm_sq = nx * nx + ny * ny + nz * nz;
418 if (n_norm_sq <= std::numeric_limits<double>::epsilon()) {
419 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
420 "Degenerated master triangle plane");
421 }
422
423 const double px = point[0] - ax;
424 const double py = point[1] - ay;
425 const double pz = point[2] - az;
426 const double signed_distance = (px * nx + py * ny + pz * nz) / n_norm_sq;
427 projected_point[0] = point[0] - signed_distance * nx;
428 projected_point[1] = point[1] - signed_distance * ny;
429 projected_point[2] = point[2] - signed_distance * nz;
430
431 MoFEMFunctionReturn(0);
432 }
433
434 static MoFEMErrorCode projectSlavePointsToMasterPlane(
435 const double *master_coords, const MatrixDouble &slave_global_coords,
436 MatrixDouble &master_global_coords) {
437 MoFEMFunctionBegin;
438
439 // The weak form is integrated on the slave surface. Each slave Gauss point
440 // is paired with its orthogonal projection on the infinite plane defined
441 // by the master triangle.
442 master_global_coords.resize(slave_global_coords.size1(), 3, false);
443 for (size_t gg = 0; gg != slave_global_coords.size1(); ++gg) {
444 CHKERR projectPointToMasterPlane(master_coords,
445 &slave_global_coords(gg, 0),
446 &master_global_coords(gg, 0));
447 }
448
449 MoFEMFunctionReturn(0);
450 }
451
452 static MoFEMErrorCode refineQuadrature(
453 FlatPrismElementForcesAndSourcesCore &fe, const int refinement_levels,
454 MatrixDouble &ref_gauss_pts) {
455 MoFEMFunctionBegin;
456
457 constexpr int num_nodes = 3;
458 moab::Core moab_ref;
459 double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
460 EntityHandle nodes[num_nodes];
461 for (int nn = 0; nn != num_nodes; ++nn)
462 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
463
464 EntityHandle tri;
465 CHKERR moab_ref.create_element(MBTRI, nodes, num_nodes, tri);
466 MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
467 MoFEM::Interface &m_field_ref = mofem_ref_core;
468
469 {
470 Range tris(tri, tri);
471 Range edges;
472 // Force MOAB to materialise the triangle edges before the bit-ref
473 // database is seeded. The working reference-refinement snippets all do
474 // this, and the side-number lookup below verifies that MOAB sees a
475 // consistent local edge numbering on the temporary reference triangle.
476 CHKERR m_field_ref.get_moab().get_adjacencies(
477 tris, 1, true, edges, moab::Interface::UNION);
478 for (auto edge : edges) {
479 int side_number = -1, sense = 0, offset = 0;
480 CHKERR moab_ref.side_number(tri, edge, side_number, sense, offset);
481 }
482 CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevel(
483 tris, BitRefLevel().set(0), false, VERBOSE);
484 }
485
486 auto *m_ref = m_field_ref.getInterface<MeshRefinement>();
487 for (int ll = 0; ll != refinement_levels; ++ll) {
488 Range tris;
489 CHKERR m_field_ref.getInterface<BitRefManager>()
490 ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
491 BitRefLevel().set(), MBTRI, tris);
492
493 Range edges;
494 CHKERR moab_ref.get_adjacencies(tris, 1, true, edges,
495 moab::Interface::UNION);
496
497 const BitRefLevel child_bit = BitRefLevel().set(ll + 1);
498 CHKERR m_ref->addVerticesInTheMiddleOfEdges(edges, child_bit);
499 CHKERR m_ref->refineTris(tris, child_bit, QUIET, false);
500 }
501
502 Range tris;
503 CHKERR m_field_ref.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
504 BitRefLevel().set(refinement_levels), BitRefLevel().set(), MBTRI, tris);
505
506 MatrixDouble ref_coords(tris.size(), 9, false);
507 int tt = 0;
508 for (auto tit = tris.begin(); tit != tris.end(); ++tit, ++tt) {
509 int num_conn;
510 const EntityHandle *conn;
511 CHKERR moab_ref.get_connectivity(*tit, conn, num_conn, false);
512 CHKERR moab_ref.get_coords(conn, num_conn, &ref_coords(tt, 0));
513 }
514
515 const size_t nb_gauss_pts = fe.gaussPts.size2();
516 MatrixDouble shape_n(nb_gauss_pts, 3, false);
517 CHKERR ShapeMBTRI(&shape_n(0, 0), &fe.gaussPts(0, 0), &fe.gaussPts(1, 0),
518 nb_gauss_pts);
519 ref_gauss_pts.resize(3, nb_gauss_pts * ref_coords.size1(), false);
520 int gg = 0;
521 for (size_t rr = 0; rr != ref_coords.size1(); ++rr) {
522 double *tri_coords = &ref_coords(rr, 0);
523 FTensor::Tensor1<double, 3> t_normal;
524 CHKERR Tools::getTriNormal(tri_coords, &t_normal(0));
525 const double det = t_normal.l2();
526 for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
527 for (int dd = 0; dd != 2; ++dd) {
528 ref_gauss_pts(dd, gg) =
529 shape_n(ggg, 0) * tri_coords[3 * 0 + dd] +
530 shape_n(ggg, 1) * tri_coords[3 * 1 + dd] +
531 shape_n(ggg, 2) * tri_coords[3 * 2 + dd];
532 }
533 ref_gauss_pts(2, gg) = fe.gaussPts(2, ggg) * det;
534 }
535 }
536
537 MoFEMFunctionReturn(0);
538 }
539
540 static inline std::map<std::pair<int, int>, MatrixDouble> mapRefCoords;
541};
542
543template <int Tensor_Dim>
544struct OpCalculateTraceVectorFieldValues
545 : public OpCalculateVectorFieldValues_General<Tensor_Dim, MatrixDouble> {
546
547 using Base = OpCalculateVectorFieldValues_General<Tensor_Dim, MatrixDouble>;
548 using Base::Base;
549
550 MoFEMErrorCode doWork(int side, EntityType type,
551 EntitiesFieldData::EntData &data) override {
552 MoFEMFunctionBeginHot;
553 if (!isOnPrismTieTrace(side, type))
554 MoFEMFunctionReturnHot(0);
555 CHKERR Base::doWork(side, type, data);
556 MoFEMFunctionReturnHot(0);
557 }
558};
559
560using TiePrismRhsBase = FormsIntegrators<
561 FlatPrismElementForcesAndSourcesCore::UserDataOperator>::Assembly<A>::OpBase;
562
563struct OpTiePrismRhs : public TiePrismRhsBase {
564
565 OpTiePrismRhs(boost::shared_ptr<MatrixDouble> disp_ptr)
566 : TiePrismRhsBase("LAMBDA", "LAMBDA", OPROW, nullptr),
567 dispPtr(disp_ptr) {
568 if (!dispPtr)
569 THROW_MESSAGE("Pointers are not set");
570 }
571
572 inline double getSlaveFaceArea() { return getAreaF3(); }
573
574protected:
575 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data) override {
576 MoFEMFunctionBeginHot;
577
578 const auto *fe =
579 dynamic_cast<const FlatPrismElementForcesAndSourcesCore *>(ptrFE);
580 if (!fe)
581 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
582 "TIE prism RHS operator is attached to a wrong finite element");
583
584 const auto nb_gauss_pts = getGaussPts().size2();
585 if (!nb_gauss_pts || nb_gauss_pts % 2) {
586 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
587 "Expected concatenated [slave|master] TIE Gauss points");
588 }
589 const auto nb_slave_gauss = nb_gauss_pts / 2;
590 if (dispPtr->size1() != nb_gauss_pts) {
591 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
592 "TIE displacement values at Gauss points were not prepared "
593 "before RHS assembly");
594 }
595
596 if (nbRows % SPACE_DIM)
597 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
598 "Unexpected number of LAMBDA dofs on TIE prism rows");
599
600 auto &row_n = row_data.getN();
601 if (!row_n.size1() || !row_n.size2())
602 MoFEMFunctionReturnHot(0);
603 if (row_n.size1() != nb_gauss_pts)
604 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
605 "Unexpected number of Gauss points on TIE multiplier block");
606
607 const int nb_row_base = nbRows / SPACE_DIM;
608 if (static_cast<size_t>(nb_row_base) > row_n.size2())
609 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
610 "More TIE multiplier bases than shape functions");
611 if (dispPtr->size2() != SPACE_DIM)
612 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
613 "Unexpected TIE displacement matrix size: size2=%d expected=%d",
614 static_cast<int>(dispPtr->size2()), SPACE_DIM);
615 if (locF.size() != static_cast<size_t>(nbRows))
616 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
617 "Unexpected TIE RHS block size: rowSide=%d rowType=%d nbRows=%d "
618 "locF.size()=%d row_n=(%d,%d) disp=(%d,%d)",
619 rowSide, rowType, nbRows, static_cast<int>(locF.size()),
620 static_cast<int>(row_n.size1()), static_cast<int>(row_n.size2()),
621 static_cast<int>(dispPtr->size1()),
622 static_cast<int>(dispPtr->size2()));
623
624 auto t_w = getFTensor0IntegrationWeight();
625 for (size_t gg = 0; gg != nb_slave_gauss; ++gg) {
626 const double slave_alpha = getSlaveFaceArea() * t_w;
627 auto t_row_base = row_data.getFTensor0N(gg, 0);
628 for (int rr = 0; rr != nb_row_base; ++rr) {
629 const double row_base = t_row_base;
630 if (row_base != 0.) {
631 auto t_row_rhs = getFTensor1FromPtr<SPACE_DIM>(&locF[SPACE_DIM * rr]);
632 for (int dd = 0; dd != SPACE_DIM; ++dd) {
633 t_row_rhs(dd) += row_base * slave_alpha *
634 ((*dispPtr)(gg, dd) -
635 (*dispPtr)(gg + nb_slave_gauss, dd));
636 }
637 }
638 ++t_row_base;
639 }
640
641 ++t_w;
642 }
643
644 MoFEMFunctionReturnHot(0);
645 }
646
647private:
648 boost::shared_ptr<MatrixDouble> dispPtr;
649};
650
651struct OpTiePrismRhsU : public TiePrismRhsBase {
652
653 OpTiePrismRhsU(const std::string &field_name,
654 boost::shared_ptr<MatrixDouble> lambda_ptr)
655 : TiePrismRhsBase(field_name, field_name, OPROW, nullptr),
656 lambdaPtr(lambda_ptr) {
657 if (!lambdaPtr)
658 THROW_MESSAGE("Pointers are not set");
659 }
660
661 inline double getSlaveFaceArea() { return getAreaF3(); }
662
663protected:
664 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data) override {
665 MoFEMFunctionBeginHot;
666
667 const auto *fe =
668 dynamic_cast<const FlatPrismElementForcesAndSourcesCore *>(ptrFE);
669 if (!fe)
670 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
671 "TIE prism RHS operator is attached to a wrong finite element");
672
673 const auto nb_gauss_pts = getGaussPts().size2();
674 if (!nb_gauss_pts || nb_gauss_pts % 2) {
675 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
676 "Expected concatenated [slave|master] TIE Gauss points");
677 }
678 const auto nb_slave_gauss = nb_gauss_pts / 2;
679 if (lambdaPtr->size1() != nb_gauss_pts) {
680 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
681 "TIE multiplier values at Gauss points were not prepared "
682 "before RHS assembly");
683 }
684
685 if (nbRows % SPACE_DIM)
686 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
687 "Unexpected number of TIE prism vector dofs on rows");
688
689 if (!isOnPrismTieTrace(rowSide, rowType))
690 MoFEMFunctionReturnHot(0);
691
692 auto &row_n = row_data.getN();
693 if (!row_n.size1() || !row_n.size2())
694 MoFEMFunctionReturnHot(0);
695 if (row_n.size1() != nb_gauss_pts)
696 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
697 "Unexpected number of Gauss points on TIE displacement block");
698
699 const int nb_row_base = nbRows / SPACE_DIM;
700 if (static_cast<size_t>(nb_row_base) > row_n.size2())
701 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
702 "More TIE displacement bases than shape functions");
703 if (lambdaPtr->size2() != SPACE_DIM)
704 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
705 "Unexpected TIE multiplier matrix size: size2=%d expected=%d",
706 static_cast<int>(lambdaPtr->size2()), SPACE_DIM);
707 if (locF.size() != static_cast<size_t>(nbRows))
708 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
709 "Unexpected TIE U RHS block size: rowSide=%d rowType=%d "
710 "nbRows=%d locF.size()=%d row_n=(%d,%d) lambda=(%d,%d)",
711 rowSide, rowType, nbRows, static_cast<int>(locF.size()),
712 static_cast<int>(row_n.size1()), static_cast<int>(row_n.size2()),
713 static_cast<int>(lambdaPtr->size1()),
714 static_cast<int>(lambdaPtr->size2()));
715
716 auto t_w = getFTensor0IntegrationWeight();
717 for (size_t gg = 0; gg != nb_slave_gauss; ++gg) {
718 const double slave_alpha = getSlaveFaceArea() * t_w;
719 for (int rr = 0; rr != nb_row_base; ++rr) {
720 const bool base_on_slave = isOnSlaveTieRowBase(rr, rowSide, rowType);
721 const double row_sign = base_on_slave ? 1.0 : -1.0;
722 const auto row_gg = base_on_slave ? gg : gg + nb_slave_gauss;
723 if (row_gg >= row_n.size1())
724 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
725 "TIE U RHS Gauss index out of range: rowSide=%d rowType=%d "
726 "rr=%d row_gg=%d row_n.size1()=%d nb_slave_gauss=%d",
727 rowSide, rowType, rr, static_cast<int>(row_gg),
728 static_cast<int>(row_n.size1()),
729 static_cast<int>(nb_slave_gauss));
730 auto t_row_base = row_data.getFTensor0N(row_gg, rr);
731 const double row_base = t_row_base;
732 if (row_base != 0.) {
733 auto t_row_rhs = getFTensor1FromPtr<SPACE_DIM>(&locF[SPACE_DIM * rr]);
734 for (int dd = 0; dd != SPACE_DIM; ++dd)
735 t_row_rhs(dd) +=
736 row_sign * row_base * slave_alpha * (*lambdaPtr)(gg, dd);
737 }
738 ++t_row_base;
739 }
740
741 ++t_w;
742 }
743
744 MoFEMFunctionReturnHot(0);
745 }
746
747private:
748 boost::shared_ptr<MatrixDouble> lambdaPtr;
749};
750
751using TiePrismLhsBase = FormsIntegrators<
752 FlatPrismElementForcesAndSourcesCore::UserDataOperator>::Assembly<A>::OpBase;
753
754struct OpTiePrismLhs : public TiePrismLhsBase {
755
756 OpTiePrismLhs(const std::string &row_field_name,
757 const std::string &col_field_name,
758 const bool assemble_transpose = false)
759 : TiePrismLhsBase(row_field_name, col_field_name, OPROWCOL, nullptr),
760 rowFieldName(row_field_name), colFieldName(col_field_name) {
761 sYmm = false;
762 this->assembleTranspose = assemble_transpose;
763 }
764
765 inline double getSlaveFaceArea() { return getAreaF3(); }
766
767protected:
768 std::string rowFieldName;
769 std::string colFieldName;
770
771 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
772 EntitiesFieldData::EntData &col_data) override {
773 MoFEMFunctionBeginHot;
774
775 const auto *fe =
776 dynamic_cast<const FlatPrismElementForcesAndSourcesCore *>(ptrFE);
777 if (!fe)
778 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
779 "TIE prism operator is attached to a wrong finite element");
780
781 if (nbRows % SPACE_DIM || nbCols % SPACE_DIM)
782 SETERRQ(
783 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
784 "Unexpected number of TIE prism vector dofs on rows or columns: "
785 "rowField=%s colField=%s rowSide=%d colSide=%d rowType=%d colType=%d "
786 "nbRows=%d nbCols=%d",
787 rowFieldName.c_str(), colFieldName.c_str(), rowSide, colSide,
788 rowType, colType, nbRows, nbCols);
789
790 if (colFieldName != "LAMBDA")
791 MoFEMFunctionReturnHot(0);
792
793 const auto nb_gauss_pts = getGaussPts().size2();
794 if (!nb_gauss_pts || nb_gauss_pts % 2) {
795 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
796 "Expected concatenated [slave|master] TIE Gauss points");
797 }
798 const auto nb_slave_gauss = nb_gauss_pts / 2;
799 const int nb_row_base = nbRows / SPACE_DIM;
800 const int nb_col_base = nbCols / SPACE_DIM;
801
802 auto &col_n = col_data.getN();
803 if (!col_n.size1() || !col_n.size2())
804 MoFEMFunctionReturnHot(0);
805
806 auto &row_n = row_data.getN();
807 if (!row_n.size1() || !row_n.size2())
808 MoFEMFunctionReturnHot(0);
809
810 if (row_n.size1() != nb_gauss_pts || col_n.size1() != nb_gauss_pts)
811 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
812 "Unexpected number of Gauss points on TIE prism blocks");
813
814 if (static_cast<size_t>(nb_row_base) > row_n.size2() ||
815 static_cast<size_t>(nb_col_base) > col_n.size2())
816 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
817 "More TIE prism bases than shape functions");
818 if (locMat.size1() != static_cast<size_t>(nbRows) ||
819 locMat.size2() != static_cast<size_t>(nbCols))
820 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
821 "Unexpected TIE LHS block size: rowField=%s colField=%s "
822 "rowSide=%d colSide=%d rowType=%d colType=%d nbRows=%d "
823 "nbCols=%d locMat=(%d,%d) row_n=(%d,%d) col_n=(%d,%d)",
824 rowFieldName.c_str(), colFieldName.c_str(), rowSide, colSide,
825 rowType, colType, nbRows, nbCols,
826 static_cast<int>(locMat.size1()), static_cast<int>(locMat.size2()),
827 static_cast<int>(row_n.size1()), static_cast<int>(row_n.size2()),
828 static_cast<int>(col_n.size1()), static_cast<int>(col_n.size2()));
829
830 if (!isOnPrismTieTrace(rowSide, rowType))
831 MoFEMFunctionReturnHot(0);
832
833 auto t_w = getFTensor0IntegrationWeight();
834 for (size_t gg = 0; gg != nb_slave_gauss; ++gg) {
835 const double slave_alpha = getSlaveFaceArea() * t_w;
836 for (int rr = 0; rr != nb_row_base; ++rr) {
837 const bool base_on_slave = isOnSlaveTieRowBase(rr, rowSide, rowType);
838 const double row_sign = base_on_slave ? 1.0 : -1.0;
839 const auto row_gg = base_on_slave ? gg : gg + nb_slave_gauss;
840 if (row_gg >= row_n.size1())
841 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
842 "TIE LHS Gauss index out of range: rowField=%s colField=%s "
843 "rowSide=%d rowType=%d rr=%d row_gg=%d row_n.size1()=%d "
844 "nb_slave_gauss=%d",
845 rowFieldName.c_str(), colFieldName.c_str(), rowSide, rowType,
846 rr, static_cast<int>(row_gg), static_cast<int>(row_n.size1()),
847 static_cast<int>(nb_slave_gauss));
848 auto t_row_base = row_data.getFTensor0N(row_gg, rr);
849 auto t_col_base = col_data.getFTensor0N(gg, 0);
850 for (int cc = 0; cc != nb_col_base; ++cc) {
851 const double value = row_sign * slave_alpha * t_row_base * t_col_base;
852 for (int dd = 0; dd != SPACE_DIM; ++dd)
853 locMat(SPACE_DIM * rr + dd, SPACE_DIM * cc + dd) += value;
854 ++t_col_base;
855 }
856 ++t_row_base;
857 }
858 ++t_w;
859 }
860
861 MoFEMFunctionReturnHot(0);
862 }
863};
864
865inline MoFEMErrorCode logTieDisplacementNorms(
866 MoFEM::Interface &mField, const boost::shared_ptr<CommonData> &tie_data,
867 const std::string &field_name) {
868 MoFEMFunctionBegin;
869
870 if (!tie_data || !tie_data->hasTieConstraints || !tie_data->tieSlaveFaces ||
871 !tie_data->tieMasterFaces) {
872 return 0;
873 }
874
875 auto &moab = mField.get_moab();
876 auto comm_interface = mField.getInterface<CommInterface>();
877 auto *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
878
879 auto closure_with_faces = [&](const Range &faces, Range &closure) {
880 MoFEMFunctionBegin;
881 closure = faces;
882 Range edges, verts;
883 CHKERR moab.get_adjacencies(faces, 1, true, edges, moab::Interface::UNION);
884 CHKERR moab.get_adjacencies(faces, 0, true, verts, moab::Interface::UNION);
885 closure.merge(edges);
886 closure.merge(verts);
887 MoFEMFunctionReturn(0);
888 };
889
890 auto calculate_norms = [&](const Range &ents, double &l2_norm, double &linf,
891 unsigned int &nb_dofs) {
892 MoFEMFunctionBegin;
893 double sum2 = 0;
894 double max_abs = 0;
895 unsigned int count = 0;
896 for (auto ent : ents) {
897 if (pcomm) {
898 unsigned char pstatus;
899 CHKERR moab.tag_get_data(pcomm->pstatus_tag(), &ent, 1, &pstatus);
900 if (pstatus & PSTATUS_NOT_OWNED)
901 continue;
902 }
903 for (_IT_GET_DOFS_FIELD_BY_NAME_AND_ENT_FOR_LOOP_(mField, field_name, ent,
904 dit)) {
905 const double value = (*dit)->getFieldData();
906 sum2 += value * value;
907 max_abs = std::max(max_abs, std::abs(value));
908 ++count;
909 }
910 }
911 double global_sum2 = 0;
912 double global_max_abs = 0;
913 unsigned long long global_count = 0;
914 const auto local_count = static_cast<unsigned long long>(count);
915 MPI_Allreduce(&sum2, &global_sum2, 1, MPI_DOUBLE, MPI_SUM,
916 mField.get_comm());
917 MPI_Allreduce(&max_abs, &global_max_abs, 1, MPI_DOUBLE, MPI_MAX,
918 mField.get_comm());
919 MPI_Allreduce(&local_count, &global_count, 1, MPI_UNSIGNED_LONG_LONG,
920 MPI_SUM, mField.get_comm());
921 l2_norm = std::sqrt(global_sum2);
922 linf = global_max_abs;
923 nb_dofs = static_cast<unsigned int>(global_count);
924 MoFEMFunctionReturn(0);
925 };
926
927 Range slave_closure, master_closure;
928 CHKERR closure_with_faces(*tie_data->tieSlaveFaces, slave_closure);
929 CHKERR closure_with_faces(*tie_data->tieMasterFaces, master_closure);
930 CHKERR comm_interface->synchroniseEntities(slave_closure);
931 CHKERR comm_interface->synchroniseEntities(master_closure);
932
933 double slave_l2 = 0, slave_linf = 0, master_l2 = 0, master_linf = 0;
934 unsigned int slave_nb_dofs = 0, master_nb_dofs = 0;
935 CHKERR calculate_norms(slave_closure, slave_l2, slave_linf, slave_nb_dofs);
936 CHKERR calculate_norms(master_closure, master_l2, master_linf, master_nb_dofs);
937
938 MOFEM_LOG("WORLD", Sev::inform)
939 << "TIE solved displacement norms on slave face closure: dofs="
940 << slave_nb_dofs << " l2=" << slave_l2 << " linf=" << slave_linf;
941 MOFEM_LOG("WORLD", Sev::inform)
942 << "TIE solved displacement norms on master face closure: dofs="
943 << master_nb_dofs << " l2=" << master_l2 << " linf=" << master_linf;
944
945 MoFEMFunctionReturn(0);
946}
947
948inline MoFEMErrorCode initializeTieConstraints(
949 MoFEM::Interface &mField, boost::shared_ptr<CommonData> &tie_data,
950 const std::string &field_name, FieldApproximationBase base, int order) {
951 MoFEMFunctionBegin;
952 auto simple = mField.getInterface<Simple>();
953 auto meshset_mng = mField.getInterface<MeshsetsManager>();
954 auto comm_interface = mField.getInterface<CommInterface>();
955
956 auto tie_master =
957 meshset_mng->getCubitMeshsetPtr(std::regex("^TIE_MASTER.*"));
958 auto tie_slave = meshset_mng->getCubitMeshsetPtr(std::regex("^TIE_SLAVE.*"));
959
960 const bool have_master_meshset = !tie_master.empty();
961 const bool have_slave_meshset = !tie_slave.empty();
962
963 if (!have_master_meshset && !have_slave_meshset) {
964 MoFEMFunctionReturnHot(0);
965 }
966
967 auto anyRank = [&](const bool local_value) {
968 int global_value = 0;
969 const int local_int = static_cast<int>(local_value);
970 MPI_Allreduce(&local_int, &global_value, 1, MPI_INT, MPI_MAX,
971 mField.get_comm());
972 return global_value != 0;
973 };
974
975 if (SPACE_DIM != 3) {
976 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
977 "TIE prism constraints are supported only in 3D");
978 }
979
980 if (!tie_data)
981 tie_data = boost::make_shared<CommonData>();
982
983 tie_data->tieRefLevel = 1;
984 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-tie_ref_level",
985 &tie_data->tieRefLevel, PETSC_NULLPTR);
986 if (tie_data->tieRefLevel < 0) {
987 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
988 "-tie_ref_level has to be >= 0");
989 }
990
991 tie_data->dispFieldName = field_name;
992 tie_data->dispFieldBase = base;
993 tie_data->dispFieldOrder = order;
994
995 if (!have_master_meshset || !have_slave_meshset) {
996 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
997 "Both TIE_MASTER* and TIE_SLAVE* meshsets are required to "
998 "initialise TIE constraints");
999 }
1000 // FIXME: this might not be necessary
1001 Range local_prisms;
1002 CHKERR mField.get_moab().get_entities_by_type(0, MBPRISM, local_prisms, true);
1003 Range local_prism_faces;
1004 if (!local_prisms.empty()) {
1005 CHKERR mField.get_moab().get_adjacencies(local_prisms, 2, false,
1006 local_prism_faces,
1007 moab::Interface::UNION);
1008 local_prism_faces = local_prism_faces.subset_by_type(MBTRI);
1009 }
1010
1011 auto collectTieFaces = [&](const auto &meshsets,
1012 boost::shared_ptr<Range> &face_range) {
1013 MoFEMFunctionBegin;
1014 face_range = boost::make_shared<Range>();
1015 for (const auto &meshset_ptr : meshsets) {
1016 Range tris;
1017 // TIE blocksets written to file can expose their surface triangles through
1018 // nested child sets after reload, so collect them recursively.
1019 CHKERR mField.get_moab().get_entities_by_type(meshset_ptr->meshset, MBTRI,
1020 tris, true);
1021 tris = intersect(tris, local_prism_faces);
1022 face_range->merge(tris);
1023 }
1024 MoFEMFunctionReturn(0);
1025 };
1026
1027 CHKERR collectTieFaces(tie_master, tie_data->tieMasterFaces);
1028 CHKERR collectTieFaces(tie_slave, tie_data->tieSlaveFaces);
1029
1030 if (tie_data->tieMasterFaces->empty() != tie_data->tieSlaveFaces->empty()) {
1031 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1032 "Local TIE_MASTER/TIE_SLAVE meshsets are inconsistent");
1033 }
1034
1035 tie_data->tiePrisms = boost::make_shared<Range>();
1036 Range master_adjacent_prisms, slave_adjacent_prisms;
1037 CHKERR mField.get_moab().get_adjacencies(*tie_data->tieSlaveFaces, 3, false,
1038 slave_adjacent_prisms,
1039 moab::Interface::UNION);
1040 CHKERR mField.get_moab().get_adjacencies(*tie_data->tieMasterFaces, 3, false,
1041 master_adjacent_prisms,
1042 moab::Interface::UNION);
1043 slave_adjacent_prisms = slave_adjacent_prisms.subset_by_type(MBPRISM);
1044 master_adjacent_prisms = master_adjacent_prisms.subset_by_type(MBPRISM);
1045
1046 if (!tie_data->tieSlaveFaces->empty() && slave_adjacent_prisms.empty()) {
1047 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1048 "No prism elements are adjacent to TIE_SLAVE meshsets");
1049 }
1050
1051 if (slave_adjacent_prisms.size() != master_adjacent_prisms.size()) {
1052 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1053 "Different numbers of prisms are adjacent to TIE_SLAVE and "
1054 "TIE_MASTER meshsets");
1055 }
1056
1057 *tie_data->tiePrisms = slave_adjacent_prisms;
1058
1059 const bool have_tie_prisms = anyRank(!tie_data->tiePrisms->empty());
1060 if (!have_tie_prisms) {
1061 MOFEM_LOG("WORLD", Sev::warning)
1062 << "TIE_MASTER/TIE_SLAVE meshsets found but no prism elements are adjacent to them. "
1063 << "TIE constraints cannot be initialised.";
1064 MoFEMFunctionReturnHot(0);
1065 }
1066
1067 if (!tie_data->tiePrisms->empty()) {
1068 MOFEM_LOG("WORLD", Sev::inform)
1069 << "Initialising TIE constraints on " << tie_data->tiePrisms->size()
1070 << " prism elements adjacent to the TIE slave/master surfaces.";
1071 }
1072
1073 if (!mField.check_field("LAMBDA"))
1074 CHKERR mField.add_field("LAMBDA", H1, base, SPACE_DIM);
1075
1076 Range tie_prism_trace_faces = *tie_data->tieSlaveFaces;
1077 tie_prism_trace_faces.merge(*tie_data->tieMasterFaces);
1078 Range tie_prism_trace_edges;
1079 CHKERR mField.get_moab().get_adjacencies(tie_prism_trace_faces, 1, false,
1080 tie_prism_trace_edges,
1081 moab::Interface::UNION);
1082 Range tie_prism_edges;
1083 CHKERR mField.get_moab().get_adjacencies(*tie_data->tiePrisms, 1, false,
1084 tie_prism_edges,
1085 moab::Interface::UNION);
1086 const auto tie_prism_lateral_edges =
1087 subtract(tie_prism_edges, tie_prism_trace_edges);
1088
1089 CHKERR mField.add_ents_to_field_by_type(*tie_data->tiePrisms, MBPRISM,
1090 field_name);
1091 CHKERR mField.add_ents_to_field_by_type(*tie_data->tiePrisms, MBPRISM,
1092 "GEOMETRY");
1093 CHKERR mField.set_field_order(tie_prism_trace_edges, "GEOMETRY", 2);
1094
1095 // LAMBDA is a classical surface multiplier. It lives only on the slave
1096 // interface, while the prism FE uses that trace space to constrain slave and
1097 // master displacement traces. Work on the already collected recursive slave
1098 // face range directly, since reloaded Cubit blocksets can expose triangles
1099 // only through nested child sets.
1100 CHKERR mField.add_ents_to_field_by_type(*tie_data->tieSlaveFaces, MBTRI,
1101 "LAMBDA");
1102 CHKERR comm_interface->synchroniseFieldEntities("LAMBDA", 0);
1103
1104 Range verts;
1105 CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, verts, true);
1106 Range faces;
1107 CHKERR mField.get_moab().get_entities_by_type(0, MBTRI, faces, true);
1108 Range edges;
1109 CHKERR mField.get_moab().get_adjacencies(faces, 1, true, edges,
1110 moab::Interface::UNION);
1111 CHKERR mField.add_ents_to_field_by_type(verts, MBVERTEX, "U");
1112 CHKERR comm_interface->synchroniseEntities(edges);
1113 CHKERR mField.add_ents_to_field_by_type(edges, MBEDGE, "U");
1114 CHKERR mField.add_ents_to_field_by_type(faces, MBTRI, "U");
1115 CHKERR comm_interface->synchroniseFieldEntities("U", 0);
1116
1117
1118 CHKERR mField.set_field_order(0, MBVERTEX, field_name, 1);
1119 CHKERR mField.set_field_order(0, MBEDGE, field_name, order);
1120 CHKERR mField.set_field_order(0, MBTRI, field_name, order);
1121 CHKERR mField.set_field_order(tie_prism_lateral_edges, field_name, 1);
1122 CHKERR mField.set_field_order(0, MBVERTEX, "LAMBDA", 1);
1123 CHKERR mField.set_field_order(0, MBEDGE, "LAMBDA", order - 1);
1124 CHKERR mField.set_field_order(0, MBTRI, "LAMBDA", order - 1);
1125
1126 CHKERR mField.build_fields();
1127
1128 if (!mField.check_finite_element(tie_data->tiePrismFeName)) {
1129 CHKERR mField.add_finite_element(tie_data->tiePrismFeName);
1130 CHKERR mField.modify_finite_element_add_field_row(tie_data->tiePrismFeName,
1131 field_name);
1132 CHKERR mField.modify_finite_element_add_field_col(tie_data->tiePrismFeName,
1133 field_name);
1134 CHKERR mField.modify_finite_element_add_field_row(tie_data->tiePrismFeName,
1135 "LAMBDA");
1136 CHKERR mField.modify_finite_element_add_field_col(tie_data->tiePrismFeName,
1137 "LAMBDA");
1138 CHKERR mField.modify_finite_element_add_field_data(tie_data->tiePrismFeName,
1139 field_name);
1140 CHKERR mField.modify_finite_element_add_field_data(tie_data->tiePrismFeName,
1141 "LAMBDA");
1142 CHKERR mField.modify_finite_element_add_field_data(tie_data->tiePrismFeName,
1143 "GEOMETRY");
1144 }
1145
1146 CHKERR mField.modify_finite_element_add_field_row(simple->getDomainFEName(),
1147 "LAMBDA");
1148 CHKERR mField.modify_finite_element_add_field_col(simple->getDomainFEName(),
1149 "LAMBDA");
1150 CHKERR mField.modify_finite_element_add_field_data(simple->getDomainFEName(),
1151 "LAMBDA");
1152
1153 CHKERR mField.add_ents_to_finite_element_by_type(*tie_data->tiePrisms,
1154 MBPRISM,
1155 tie_data->tiePrismFeName);
1156 CHKERR mField.build_finite_elements(simple->getDomainFEName());
1157 CHKERR mField.build_finite_elements(tie_data->tiePrismFeName);
1158
1159 auto &other_fes = simple->getOtherFiniteElements();
1160 other_fes.push_back(tie_data->tiePrismFeName);
1161
1162 CHKERR DMMoFEMAddElement(simple->getDM(), tie_data->tiePrismFeName);
1163 CHKERR simple->reSetUp(true);
1164
1165 tie_data->hasTieConstraints = true;
1166 MoFEMFunctionReturn(0);
1167}
1168
1169inline MoFEMErrorCode setupTieConstraintOperators(
1170 MoFEM::Interface &mField, boost::shared_ptr<CommonData> tie_data) {
1171 MoFEMFunctionBegin;
1172
1173 if (!tie_data || !tie_data->hasTieConstraints) {
1174 MoFEMFunctionReturnHot(0);
1175 }
1176
1177 if (!tie_data->tiePrismRhsFE) {
1178 tie_data->tiePrismRhsFE =
1179 boost::make_shared<FlatPrismElementForcesAndSourcesCore>(mField);
1180 tie_data->tiePrismRhsFE->getRuleHook = [](int, int, int) { return -1; };
1181 tie_data->tiePrismRhsFE->setRuleHook =
1182 SetIntegrationOnTiePrismFaces(tie_data->tieRefLevel);
1183 tie_data->tiePrismRhsFE->getOpPtrVector().push_back(
1184 new OpCalculateTraceVectorFieldValues<SPACE_DIM>(
1185 "LAMBDA", tie_data->lambdaAtGaussPts, MBVERTEX));
1186 tie_data->tiePrismRhsFE->getOpPtrVector().push_back(
1187 new OpTiePrismRhsU(tie_data->dispFieldName,
1188 tie_data->lambdaAtGaussPts));
1189 tie_data->tiePrismRhsFE->getOpPtrVector().push_back(
1190 new OpCalculateTraceVectorFieldValues<SPACE_DIM>(
1191 tie_data->dispFieldName, tie_data->dispAtGaussPts, MBVERTEX));
1192 tie_data->tiePrismRhsFE->getOpPtrVector().push_back(
1193 new OpTiePrismRhs(tie_data->dispAtGaussPts));
1194 }
1195
1196 if (!tie_data->tiePrismLhsFE) {
1197 tie_data->tiePrismLhsFE =
1198 boost::make_shared<FlatPrismElementForcesAndSourcesCore>(mField);
1199 tie_data->tiePrismLhsFE->getRuleHook = [](int, int, int) { return -1; };
1200 tie_data->tiePrismLhsFE->setRuleHook =
1201 SetIntegrationOnTiePrismFaces(tie_data->tieRefLevel);
1202 tie_data->tiePrismLhsFE->getOpPtrVector().push_back(
1203 new OpTiePrismLhs(tie_data->dispFieldName, "LAMBDA", true));
1204 }
1205
1206 MoFEMFunctionReturn(0);
1207}
1208
1209inline MoFEMErrorCode setupTieSolver(MoFEM::Interface &mField,
1210 boost::shared_ptr<CommonData> tie_data,
1211 SmartPetscObj<KSP> solver) {
1212 MoFEMFunctionBegin;
1213
1214 if (!tie_data || !tie_data->hasTieConstraints)
1215 MoFEMFunctionReturnHot(0);
1216
1217 CHKERR setupTieConstraintOperators(mField, tie_data);
1218
1219 DM dm;
1220 CHKERR KSPGetDM(solver, &dm);
1221 boost::shared_ptr<FEMethod> null;
1222 CHKERR DMMoFEMKSPSetComputeOperators(dm, tie_data->tiePrismFeName,
1223 tie_data->tiePrismLhsFE, null, null);
1224 CHKERR DMMoFEMKSPSetComputeRHS(dm, tie_data->tiePrismFeName,
1225 tie_data->tiePrismRhsFE, null, null);
1226
1227 MoFEMFunctionReturn(0);
1228}
1229
1230inline MoFEMErrorCode setupTieSolver(MoFEM::Interface &mField,
1231 boost::shared_ptr<CommonData> tie_data,
1232 SmartPetscObj<TS> solver) {
1233 MoFEMFunctionBegin;
1234
1235 if (!tie_data || !tie_data->hasTieConstraints)
1236 MoFEMFunctionReturnHot(0);
1237
1238 CHKERR setupTieConstraintOperators(mField, tie_data);
1239
1240 DM dm;
1241 CHKERR TSGetDM(solver, &dm);
1242 boost::shared_ptr<FEMethod> null;
1243 CHKERR DMMoFEMTSSetIJacobian(dm, tie_data->tiePrismFeName,
1244 tie_data->tiePrismLhsFE, null, null);
1245 CHKERR DMMoFEMTSSetIFunction(dm, tie_data->tiePrismFeName,
1246 tie_data->tiePrismRhsFE, null, null);
1247
1248 MoFEMFunctionReturn(0);
1249}
1250
1251} // namespace ElasticTie
1252
1253/** @brief ElasticExample extension providing reusable TIE mesh behaviour. */
1254struct ElasticTieMeshExample : public ElasticExample {
1255
1256 using ElasticExample::ElasticExample;
1257
1258protected:
1259 boost::shared_ptr<ElasticTie::CommonData> tieData = nullptr;
1260 Range physicalBoundaryFaces;
1261
1262 MoFEMErrorCode readMesh() override;
1263 MoFEMErrorCode setupProblem() override;
1264 MoFEMErrorCode setUpSolver(SmartPetscObj<KSP> solver) override;
1265 MoFEMErrorCode kspSetUpAndSolve(SmartPetscObj<KSP> solver) override;
1266 MoFEMErrorCode checkResults() override;
1267};
1268
1269inline MoFEMErrorCode ElasticTieMeshExample::readMesh() {
1270 MoFEMFunctionBegin;
1271 CHKERR ElasticTie::loadTieMesh<SPACE_DIM>(mField, physicalBoundaryFaces);
1272 MoFEMFunctionReturn(0);
1273}
1274
1275inline MoFEMErrorCode ElasticTieMeshExample::setupProblem() {
1276 MoFEMFunctionBegin;
1277 auto simple = mField.getInterface<Simple>();
1278
1279 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
1280 const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
1281 PetscInt choice_base_value = AINSWORTH;
1282 CHKERR PetscOptionsGetEList(PETSC_NULLPTR, PETSC_NULLPTR, "-base", list_bases,
1283 LASBASETOPT, &choice_base_value, PETSC_NULLPTR);
1284
1285 FieldApproximationBase base;
1286 switch (choice_base_value) {
1287 case AINSWORTH:
1288 base = AINSWORTH_LEGENDRE_BASE;
1289 break;
1290 case DEMKOWICZ:
1291 base = DEMKOWICZ_JACOBI_BASE;
1292 break;
1293 default:
1294 base = LASTBASE;
1295 break;
1296 }
1297
1298 int order = 3;
1299 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order,
1300 PETSC_NULLPTR);
1301
1302 if (SPACE_DIM == 3) {
1303 Range domain_ents;
1304 CHKERR mField.get_moab().get_entities_by_handle(simple->getMeshset(),
1305 domain_ents, true);
1306 auto tet_hex_ents = ElasticTie::getPhysicalVolumeEntities(domain_ents);
1307
1308 if (!tet_hex_ents.empty()) {
1309 EntityHandle domain_meshset;
1310 CHKERR ElasticTie::makeClosureMeshset(
1311 mField.get_moab(), tet_hex_ents, SPACE_DIM, domain_meshset);
1312 simple->getMeshset() = domain_meshset;
1313
1314 Range boundary_faces;
1315 if (mField.get_comm_size() > 1) {
1316 Skinner skin(&mField.get_moab());
1317 Range local_skin;
1318 CHKERR skin.find_skin(0, tet_hex_ents, false, local_skin);
1319 boundary_faces = intersect(physicalBoundaryFaces, local_skin);
1320 } else if (!physicalBoundaryFaces.empty()) {
1321 boundary_faces = physicalBoundaryFaces;
1322 } else {
1323 Skinner skin(&mField.get_moab());
1324 CHKERR skin.find_skin(0, tet_hex_ents, false, boundary_faces);
1325 }
1326 EntityHandle boundary_meshset;
1327 CHKERR ElasticTie::makeClosureMeshset(mField.get_moab(), boundary_faces,
1328 SPACE_DIM - 1, boundary_meshset);
1329 simple->getBoundaryMeshSet() = boundary_meshset;
1330 }
1331 }
1332
1333 CHKERR ElasticExample::setupProblem();
1334 CHKERR ElasticTie::initializeTieConstraints(mField, tieData, "U", base,
1335 order);
1336
1337 MoFEMFunctionReturn(0);
1338}
1339
1340inline MoFEMErrorCode
1346
1347inline MoFEMErrorCode
1354
1359
1360#endif // __ELASTIC_TIE_HPP__
ElasticExample.hpp
Implementation of elastic example class.
type
std::string type
Definition JsonConfigManager.cpp:28
FTENSOR_INDEX
#define FTENSOR_INDEX(DIM, I)
Definition Templates.hpp:2561
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
SPACE_DIM
constexpr int SPACE_DIM
Definition child_and_parent.cpp:17
QUIET
@ QUIET
Definition definitions.h:221
VERBOSE
@ VERBOSE
Definition definitions.h:222
FieldApproximationBase
FieldApproximationBase
approximation base
Definition definitions.h:58
LASTBASE
@ LASTBASE
Definition definitions.h:69
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
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:463
H1
@ H1
continuous field
Definition definitions.h:85
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition definitions.h:228
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
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
THROW_MESSAGE
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition definitions.h:577
order
constexpr int order
Definition dg_projection.cpp:38
ShapeMBTRI
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition fem_tools.c:182
eta
double eta
Definition free_surface.cpp:192
_IT_GET_DOFS_FIELD_BY_NAME_AND_ENT_FOR_LOOP_
#define _IT_GET_DOFS_FIELD_BY_NAME_AND_ENT_FOR_LOOP_(MFIELD, NAME, ENT, IT)
loop over all dofs from a moFEM field and particular field
Definition Interface.hpp:2023
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_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::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
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.
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
tol
double tol
Definition mesh_smoothing.cpp:23
ElasticTie::isInsideReferenceTriangle
bool isInsideReferenceTriangle(const double xi, const double eta, const double tol=1e-8)
Definition ElasticTie.hpp:197
ElasticTie::MASTER_FACE_SIDE
constexpr int MASTER_FACE_SIDE
Definition ElasticTie.hpp:23
ElasticTie::loadTieMesh
MoFEMErrorCode loadTieMesh(MoFEM::Interface &mField, Range &physicalBoundaryFaces)
Load a mesh while preserving the shared entities needed by TIE constraints.
Definition ElasticTie.hpp:89
ElasticTie::getPhysicalVolumeEntities
Range getPhysicalVolumeEntities(const Range &ents)
Definition ElasticTie.hpp:61
ElasticTie::isOnPrismTieTrace
bool isOnPrismTieTrace(const int side, const EntityType type)
Definition ElasticTie.hpp:218
ElasticTie::isOnSlaveTieRowBase
bool isOnSlaveTieRowBase(const int rr, const int row_side, const EntityType row_type)
Definition ElasticTie.hpp:223
ElasticTie::initializeTieConstraints
MoFEMErrorCode initializeTieConstraints(MoFEM::Interface &mField, boost::shared_ptr< CommonData > &tie_data, const std::string &field_name, FieldApproximationBase base, int order)
Definition ElasticTie.hpp:948
ElasticTie::isOnPrismTriFace
bool isOnPrismTriFace(const int side, const EntityType type, const int prism_face_side)
Definition ElasticTie.hpp:202
ElasticTie::getPrismTriangleFaces
MoFEMErrorCode getPrismTriangleFaces(moab::Interface &moab, const Range &ents, Range &triangle_faces)
Definition ElasticTie.hpp:67
ElasticTie::setupTieConstraintOperators
MoFEMErrorCode setupTieConstraintOperators(MoFEM::Interface &mField, boost::shared_ptr< CommonData > tie_data)
Definition ElasticTie.hpp:1169
ElasticTie::SLAVE_FACE_SIDE
constexpr int SLAVE_FACE_SIDE
Definition ElasticTie.hpp:22
ElasticTie::logTieDisplacementNorms
MoFEMErrorCode logTieDisplacementNorms(MoFEM::Interface &mField, const boost::shared_ptr< CommonData > &tie_data, const std::string &field_name)
Definition ElasticTie.hpp:865
ElasticTie::makeClosureMeshset
MoFEMErrorCode makeClosureMeshset(moab::Interface &moab, const Range &ents, const int max_dim, EntityHandle &meshset)
Definition ElasticTie.hpp:45
ElasticTie::setupTieSolver
MoFEMErrorCode setupTieSolver(MoFEM::Interface &mField, boost::shared_ptr< CommonData > tie_data, SmartPetscObj< KSP > solver)
Definition ElasticTie.hpp:1209
field_name
constexpr auto field_name
Definition poisson_2d_homogeneous.cpp:13
QUAD_2D_TABLE_SIZE
#define QUAD_2D_TABLE_SIZE
Definition quad.h:174
QUAD_2D_TABLE
static QUAD *const QUAD_2D_TABLE[]
Definition quad.h:175
ElasticExample::setupProblem
virtual MoFEMErrorCode setupProblem()
[Read mesh]
Definition ElasticExample.hpp:75
ElasticExample::kspSetUpAndSolve
virtual MoFEMErrorCode kspSetUpAndSolve(SmartPetscObj< KSP > solver)
Definition ElasticExample.hpp:212
ElasticExample::mField
MoFEM::Interface & mField
Definition ElasticExample.hpp:29
ElasticExample::ElasticExample
ElasticExample(MoFEM::Interface &m_field)
Definition ElasticExample.hpp:20
ElasticTieMeshExample
ElasticExample extension providing reusable TIE mesh behaviour.
Definition ElasticTie.hpp:1254
ElasticTieMeshExample::readMesh
MoFEMErrorCode readMesh() override
[Run problem]
Definition ElasticTie.hpp:1269
ElasticTieMeshExample::setupProblem
MoFEMErrorCode setupProblem() override
[Read mesh]
Definition ElasticTie.hpp:1275
ElasticTieMeshExample::tieData
boost::shared_ptr< ElasticTie::CommonData > tieData
Definition ElasticTie.hpp:1259
ElasticTieMeshExample::checkResults
MoFEMErrorCode checkResults() override
[Postprocess results]
Definition ElasticTie.hpp:1355
ElasticTieMeshExample::setUpSolver
MoFEMErrorCode setUpSolver(SmartPetscObj< KSP > solver) override
[Push operators to pipeline]
Definition ElasticTie.hpp:1341
ElasticTieMeshExample::kspSetUpAndSolve
MoFEMErrorCode kspSetUpAndSolve(SmartPetscObj< KSP > solver) override
Definition ElasticTie.hpp:1348
ElasticTie::CommonData::dispAtGaussPts
boost::shared_ptr< MatrixDouble > dispAtGaussPts
Definition ElasticTie.hpp:35
ElasticTie::CommonData::tiePrisms
boost::shared_ptr< Range > tiePrisms
Definition ElasticTie.hpp:34
ElasticTie::CommonData::tieSlaveFaces
boost::shared_ptr< Range > tieSlaveFaces
Definition ElasticTie.hpp:33
ElasticTie::CommonData::tiePrismLhsFE
boost::shared_ptr< FlatPrismElementForcesAndSourcesCore > tiePrismLhsFE
Definition ElasticTie.hpp:41
ElasticTie::CommonData::tiePrismFeName
std::string tiePrismFeName
Definition ElasticTie.hpp:27
ElasticTie::CommonData::lambdaAtGaussPts
boost::shared_ptr< MatrixDouble > lambdaAtGaussPts
Definition ElasticTie.hpp:37
ElasticTie::CommonData::dispFieldBase
FieldApproximationBase dispFieldBase
Definition ElasticTie.hpp:29
ElasticTie::CommonData::dispFieldName
std::string dispFieldName
Definition ElasticTie.hpp:28
ElasticTie::CommonData::tiePrismRhsFE
boost::shared_ptr< FlatPrismElementForcesAndSourcesCore > tiePrismRhsFE
Definition ElasticTie.hpp:39
ElasticTie::CommonData::tieMasterFaces
boost::shared_ptr< Range > tieMasterFaces
Definition ElasticTie.hpp:32
ElasticTie::OpCalculateTraceVectorFieldValues::Base
OpCalculateVectorFieldValues_General< Tensor_Dim, MatrixDouble > Base
Definition ElasticTie.hpp:547
ElasticTie::OpCalculateTraceVectorFieldValues::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data) override
Definition ElasticTie.hpp:550
ElasticTie::OpTiePrismLhs::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data) override
Definition ElasticTie.hpp:771
ElasticTie::OpTiePrismLhs::OpTiePrismLhs
OpTiePrismLhs(const std::string &row_field_name, const std::string &col_field_name, const bool assemble_transpose=false)
Definition ElasticTie.hpp:756
ElasticTie::OpTiePrismRhsU::OpTiePrismRhsU
OpTiePrismRhsU(const std::string &field_name, boost::shared_ptr< MatrixDouble > lambda_ptr)
Definition ElasticTie.hpp:653
ElasticTie::OpTiePrismRhsU::lambdaPtr
boost::shared_ptr< MatrixDouble > lambdaPtr
Definition ElasticTie.hpp:748
ElasticTie::OpTiePrismRhsU::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data) override
Definition ElasticTie.hpp:664
ElasticTie::OpTiePrismRhs::OpTiePrismRhs
OpTiePrismRhs(boost::shared_ptr< MatrixDouble > disp_ptr)
Definition ElasticTie.hpp:565
ElasticTie::OpTiePrismRhs::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data) override
Definition ElasticTie.hpp:575
ElasticTie::OpTiePrismRhs::dispPtr
boost::shared_ptr< MatrixDouble > dispPtr
Definition ElasticTie.hpp:648
ElasticTie::SetIntegrationOnTiePrismFaces::setBaseQuadrature
static MoFEMErrorCode setBaseQuadrature(FlatPrismElementForcesAndSourcesCore &fe, const int rule)
Definition ElasticTie.hpp:376
ElasticTie::SetIntegrationOnTiePrismFaces::SetIntegrationOnTiePrismFaces
SetIntegrationOnTiePrismFaces(int tie_ref_level=1)
Definition ElasticTie.hpp:245
ElasticTie::SetIntegrationOnTiePrismFaces::projectSlavePointsToMasterPlane
static MoFEMErrorCode projectSlavePointsToMasterPlane(const double *master_coords, const MatrixDouble &slave_global_coords, MatrixDouble &master_global_coords)
Definition ElasticTie.hpp:434
ElasticTie::SetIntegrationOnTiePrismFaces::mapRefCoords
static std::map< std::pair< int, int >, MatrixDouble > mapRefCoords
Definition ElasticTie.hpp:540
ElasticTie::SetIntegrationOnTiePrismFaces::refineQuadrature
static MoFEMErrorCode refineQuadrature(FlatPrismElementForcesAndSourcesCore &fe, const int refinement_levels, MatrixDouble &ref_gauss_pts)
Definition ElasticTie.hpp:452
ElasticTie::SetIntegrationOnTiePrismFaces::SetIntegrationOnTiePrismFaces
SetIntegrationOnTiePrismFaces(FunRule fun_rule, int tie_ref_level=1)
Definition ElasticTie.hpp:248
ElasticTie::SetIntegrationOnTiePrismFaces::projectPointToMasterPlane
static MoFEMErrorCode projectPointToMasterPlane(const double *master_coords, const double *point, double *projected_point)
Definition ElasticTie.hpp:390
ElasticTie::SetIntegrationOnTiePrismFaces::operator()
MoFEMErrorCode operator()(ForcesAndSourcesCore *fe_raw_ptr, int order_row, int order_col, int order_data)
Definition ElasticTie.hpp:251
MoFEM::CoreInterface::get_comm_size
virtual int get_comm_size() const =0
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::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
virtual MPI_Comm & get_comm() const =0
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93
QUAD_::npoints
int npoints
Definition quad.h:29
Generated by Doxygen 1.9.8 and hosted at University of Glasgow