v0.13.2
Loading...
Searching...
No Matches
VolumeElementForcesAndSourcesCoreOnContactPrismSide.cpp
Go to the documentation of this file.
1/** \file VolumeElementForcesAndSourcesCoreOnContactPrismSide.cpp
2
3\brief Implementation of volume element on side
4
5*/
6
7
8
9namespace MoFEM {
10
12 return -1;
13};
14
18
19 if (!sidePtrFE)
20 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
21 "Side volume element not set");
22
23 const EntityHandle volume_entity =
25
26 Range adj_faces;
27 CHKERR sidePtrFE->mField.getInterface<BitRefManager>()->getAdjacenciesAny(
28 volume_entity, 2, adj_faces);
29 SideNumber_multiIndex &side_table = const_cast<SideNumber_multiIndex &>(
30 numeredEntFiniteElementPtr->getSideNumberTable());
31
32 SideNumber_multiIndex &side_volume_table =
33 const_cast<SideNumber_multiIndex &>(
34 sidePtrFE->numeredEntFiniteElementPtr->getSideNumberTable());
35
36 adj_faces = adj_faces.subset_by_type(MBTRI);
37 bool face_common = false;
38 SideNumber_multiIndex::nth_index<0>::type::iterator sit;
39 SideNumber_multiIndex::nth_index<0>::type::iterator sit_volume;
40 EntityHandle face_ent;
41 for (Range::iterator it_tris = adj_faces.begin(); it_tris != adj_faces.end();
42 it_tris++) {
43 sit = side_table.get<0>().find(*it_tris);
44 sit_volume = side_volume_table.get<0>().find(*it_tris);
45 if (sit != side_table.get<0>().end()) {
46 face_ent = *it_tris;
47 face_common = true;
48 break;
49 }
50 }
51 if (!face_common) {
52 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Error no common tets");
53 }
54
55 auto contact_prism_ptr_fe =
57
58 faceSense = (*sit)->sense;
59 faceSideNumber = (*sit)->side_number;
60 int side_of_vol_number = (*sit_volume)->side_number;
61 const EntityHandle *ent_on_vol;
62 int num_nodes = 3;
63 CHKERR contact_prism_ptr_fe->mField.get_moab().get_connectivity(
64 face_ent, ent_on_vol, num_nodes, true);
65
66 fill(tetConnMap.begin(), tetConnMap.end(), -1);
67 for (int nn = 0; nn != 3; ++nn) {
68 faceConnMap[nn] = std::distance(conn, find(conn, &conn[4], ent_on_vol[nn]));
69
70 tetConnMap[faceConnMap[nn]] = nn;
71 if (faceConnMap[nn] > 3) {
72 cerr << " faceConnMap[nn] " << faceConnMap[nn] << "\n";
73 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
74 "No common node on face and element can not be found");
75 }
76 }
77
78 oppositeNode = std::distance(tetConnMap.begin(),
79 find(tetConnMap.begin(), tetConnMap.end(), -1));
80
81 const int nb_gauss_pts =
82 contact_prism_ptr_fe->getGaussPtsMasterFromEleSide().size2();
83 gaussPts.resize(4, nb_gauss_pts, false);
84 gaussPts.clear();
85
86 const EntityType tri_type = MBTRI;
87 boost::shared_ptr<EntitiesFieldData> dataH1_on_face;
88
89 if (side_of_vol_number == 3) {
90 dataH1_on_face = boost::shared_ptr<EntitiesFieldData>(
91 contact_prism_ptr_fe->getDataOnMasterFromEleSide()[H1]);
92 } else if (side_of_vol_number == 4) {
93 dataH1_on_face = boost::shared_ptr<EntitiesFieldData>(
94 contact_prism_ptr_fe->getDataOnSlaveFromEleSide()[H1]);
95 } else {
96 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
97 "Side element for contact: Wrong face for contact!");
98 }
99
100 const MatrixDouble &face_shape_funtions =
101 dataH1_on_face->dataOnEntities[MBVERTEX][0].getN(NOBASE);
102
103 const double tet_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
104 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
105 gaussPts(0, gg) =
106 face_shape_funtions(gg, 0) * tet_coords[3 * faceConnMap[0] + 0] +
107 face_shape_funtions(gg, 1) * tet_coords[3 * faceConnMap[1] + 0] +
108 face_shape_funtions(gg, 2) * tet_coords[3 * faceConnMap[2] + 0];
109 gaussPts(1, gg) =
110 face_shape_funtions(gg, 0) * tet_coords[3 * faceConnMap[0] + 1] +
111 face_shape_funtions(gg, 1) * tet_coords[3 * faceConnMap[1] + 1] +
112 face_shape_funtions(gg, 2) * tet_coords[3 * faceConnMap[2] + 1];
113 gaussPts(2, gg) =
114 face_shape_funtions(gg, 0) * tet_coords[3 * faceConnMap[0] + 2] +
115 face_shape_funtions(gg, 1) * tet_coords[3 * faceConnMap[1] + 2] +
116 face_shape_funtions(gg, 2) * tet_coords[3 * faceConnMap[2] + 2];
117 if (side_of_vol_number == 3) {
118 gaussPts(3, gg) =
119 contact_prism_ptr_fe->getGaussPtsMasterFromEleSide()(2, gg);
120 } else {
121 gaussPts(3, gg) =
122 contact_prism_ptr_fe->getGaussPtsSlaveFromEleSide()(2, gg);
123 }
124 }
126}
127
128} // namespace MoFEM
@ NOBASE
Definition: definitions.h:59
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
multi_index_container< boost::shared_ptr< SideNumber >, indexed_by< ordered_unique< member< SideNumber, EntityHandle, &SideNumber::ent > >, ordered_non_unique< composite_key< SideNumber, const_mem_fun< SideNumber, EntityType, &SideNumber::getEntType >, member< SideNumber, signed char, &SideNumber::side_number > > > > > SideNumber_multiIndex
SideNumber_multiIndex for SideNumber.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
Managing BitRefLevels.
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
ForcesAndSourcesCore * sidePtrFE
Element to integrate on the sides.
MatrixDouble gaussPts
Matrix of integration points.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.