v0.10.0
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 /* This file is part of MoFEM.
8  * MoFEM is free software: you can redistribute it and/or modify it under
9  * the terms of the GNU Lesser General Public License as published by the
10  * Free Software Foundation, either version 3 of the License, or (at your
11  * option) any later version.
12  *
13  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  * License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20 
21 namespace MoFEM {
22 
24  return -1;
25 };
26 
30 
31  if (!sidePtrFE)
32  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
33  "Side volume element not set");
34 
35  const EntityHandle volume_entity =
37 
38  Range adj_faces;
39  CHKERR sidePtrFE->mField.getInterface<BitRefManager>()->getAdjacenciesAny(
40  volume_entity, 2, adj_faces);
41  SideNumber_multiIndex &side_table = const_cast<SideNumber_multiIndex &>(
42  numeredEntFiniteElementPtr->getSideNumberTable());
43 
44  SideNumber_multiIndex &side_volume_table =
45  const_cast<SideNumber_multiIndex &>(
46  sidePtrFE->numeredEntFiniteElementPtr->getSideNumberTable());
47 
48  adj_faces = adj_faces.subset_by_type(MBTRI);
49  bool face_common = false;
50  SideNumber_multiIndex::nth_index<0>::type::iterator sit;
51  SideNumber_multiIndex::nth_index<0>::type::iterator sit_volume;
52  EntityHandle face_ent;
53  for (Range::iterator it_tris = adj_faces.begin(); it_tris != adj_faces.end();
54  it_tris++) {
55  sit = side_table.get<0>().find(*it_tris);
56  sit_volume = side_volume_table.get<0>().find(*it_tris);
57  if (sit != side_table.get<0>().end()) {
58  face_ent = *it_tris;
59  face_common = true;
60  break;
61  }
62  }
63  if (!face_common) {
64  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Error no common tets");
65  }
66 
67  auto contact_prism_ptr_fe =
68  static_cast<ContactPrismElementForcesAndSourcesCore *>(sidePtrFE);
69 
70  faceSense = (*sit)->sense;
71  faceSideNumber = (*sit)->side_number;
72  int side_of_vol_number = (*sit_volume)->side_number;
73  const EntityHandle *ent_on_vol;
74  int num_nodes = 3;
75  CHKERR contact_prism_ptr_fe->mField.get_moab().get_connectivity(
76  face_ent, ent_on_vol, num_nodes, true);
77 
78  fill(tetConnMap.begin(), tetConnMap.end(), -1);
79  for (int nn = 0; nn != 3; ++nn) {
80  faceConnMap[nn] = std::distance(conn, find(conn, &conn[4], ent_on_vol[nn]));
81 
82  tetConnMap[faceConnMap[nn]] = nn;
83  if (faceConnMap[nn] > 3) {
84  cerr << " faceConnMap[nn] " << faceConnMap[nn] << "\n";
85  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
86  "No common node on face and element can not be found");
87  }
88  }
89 
90  oppositeNode = std::distance(tetConnMap.begin(),
91  find(tetConnMap.begin(), tetConnMap.end(), -1));
92 
93  const int nb_gauss_pts =
94  contact_prism_ptr_fe->getGaussPtsMasterFromEleSide().size2();
95  gaussPts.resize(4, nb_gauss_pts, false);
96  gaussPts.clear();
97 
98  const EntityType tri_type = MBTRI;
99  boost::shared_ptr<DataForcesAndSourcesCore> dataH1_on_face;
100 
101  if (side_of_vol_number == 3) {
102  dataH1_on_face = boost::shared_ptr<DataForcesAndSourcesCore>(
103  contact_prism_ptr_fe->getDataOnMasterFromEleSide()[H1]);
104  } else if (side_of_vol_number == 4) {
105  dataH1_on_face = boost::shared_ptr<DataForcesAndSourcesCore>(
106  contact_prism_ptr_fe->getDataOnSlaveFromEleSide()[H1]);
107  } else {
108  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
109  "Side element for contact: Wrong face for contact!");
110  }
111 
112  const MatrixDouble &face_shape_funtions =
113  dataH1_on_face->dataOnEntities[MBVERTEX][0].getN(NOBASE);
114 
115  const double tet_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
116  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
117  gaussPts(0, gg) =
118  face_shape_funtions(gg, 0) * tet_coords[3 * faceConnMap[0] + 0] +
119  face_shape_funtions(gg, 1) * tet_coords[3 * faceConnMap[1] + 0] +
120  face_shape_funtions(gg, 2) * tet_coords[3 * faceConnMap[2] + 0];
121  gaussPts(1, gg) =
122  face_shape_funtions(gg, 0) * tet_coords[3 * faceConnMap[0] + 1] +
123  face_shape_funtions(gg, 1) * tet_coords[3 * faceConnMap[1] + 1] +
124  face_shape_funtions(gg, 2) * tet_coords[3 * faceConnMap[2] + 1];
125  gaussPts(2, gg) =
126  face_shape_funtions(gg, 0) * tet_coords[3 * faceConnMap[0] + 2] +
127  face_shape_funtions(gg, 1) * tet_coords[3 * faceConnMap[1] + 2] +
128  face_shape_funtions(gg, 2) * tet_coords[3 * faceConnMap[2] + 2];
129  if (side_of_vol_number == 3) {
130  gaussPts(3, gg) =
131  contact_prism_ptr_fe->getGaussPtsMasterFromEleSide()(2, gg);
132  } else {
133  gaussPts(3, gg) =
134  contact_prism_ptr_fe->getGaussPtsSlaveFromEleSide()(2, gg);
135  }
136  }
138 }
139 
140 } // namespace MoFEM
MatrixDouble gaussPts
Matrix of integration points.
ForcesAndSourcesCore * sidePtrFE
Element to integrate on the sides.
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, char, &SideNumber::side_number > > > > > SideNumber_multiIndex
SideNumber_multiIndex for SideNumber.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
constexpr int order
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
Managing BitRefLevels.
#define CHKERR
Inline error check.
Definition: definitions.h:604
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
continuous field
Definition: definitions.h:177