v0.14.0
PrismsFromSurfaceInterface.cpp
Go to the documentation of this file.
1 /** \file PrismsFromSurfaceInterface.cpp
2  * \brief Interface for creating prisms from surface elements
3  */
4 
5 
6 
7 namespace MoFEM {
8 
9 MoFEMErrorCode PrismsFromSurfaceInterface::query_interface(
10  boost::typeindex::type_index type_index, UnknownInterface **iface) const {
11  *iface = const_cast<PrismsFromSurfaceInterface *>(this);
12  return 0;
13 }
14 
15 DEPRECATED MoFEMErrorCode PrismsFromSurfaceInterface::createPrisms(
16  const Range &ents, Range &prisms, int verb) {
17  MoFEMFunctionBegin;
18  CHKERR createPrisms(ents, NO_SWAP, prisms, verb);
19  MoFEMFunctionReturn(0);
20 }
21 
22 MoFEMErrorCode PrismsFromSurfaceInterface::createPrisms(
23  const Range &ents, const SwapType swap_type, Range &prisms, int verb) {
24  MoFEMFunctionBegin;
25  Interface &m_field = cOre;
26  Range tris = ents.subset_by_type(MBTRI);
27  for (Range::iterator tit = tris.begin(); tit != tris.end(); tit++) {
28  const EntityHandle *conn;
29  int number_nodes = 0;
30  CHKERR m_field.get_moab().get_connectivity(*tit, conn, number_nodes, false);
31  double coords[3 * number_nodes];
32  CHKERR m_field.get_moab().get_coords(conn, number_nodes, coords);
33  EntityHandle prism_nodes[6];
34  for (int nn = 0; nn < 3; nn++) {
35  prism_nodes[nn] = conn[nn];
36  if (createdVertices.find(conn[nn]) != createdVertices.end()) {
37  prism_nodes[3 + nn] = createdVertices[prism_nodes[nn]];
38  } else {
39  CHKERR m_field.get_moab().create_vertex(&coords[3 * nn],
40  prism_nodes[3 + nn]);
41  createdVertices[conn[nn]] = prism_nodes[3 + nn];
42  CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(),
43  &prism_nodes[3 + nn], 1,
44  &prism_nodes[nn]);
45  }
46  }
47  switch (swap_type) {
48  case SWAP_TRI_NODE_ORDER:
49  std::swap(prism_nodes[1], prism_nodes[2]);
50  std::swap(prism_nodes[4], prism_nodes[5]);
51  break;
52  case SWAP_TOP_AND_BOT_TRI:
53  std::swap(prism_nodes[0], prism_nodes[3]);
54  std::swap(prism_nodes[1], prism_nodes[4]);
55  std::swap(prism_nodes[2], prism_nodes[5]);
56  break;
57  case NO_SWAP:
58  default:
59  break;
60  }
61  EntityHandle prism;
62  CHKERR m_field.get_moab().create_element(MBPRISM, prism_nodes, 6, prism);
63  Range edges;
64  CHKERR m_field.get_moab().get_adjacencies(&prism, 1, 1, true, edges,
65  moab::Interface::UNION);
66  Range faces;
67  CHKERR m_field.get_moab().get_adjacencies(&prism, 1, 2, true, faces,
68  moab::Interface::UNION);
69  prisms.insert(prism);
70  for (int ee = 0; ee <= 2; ee++) {
71  EntityHandle e1;
72  CHKERR m_field.get_moab().side_element(prism, 1, ee, e1);
73  EntityHandle e2;
74  CHKERR m_field.get_moab().side_element(prism, 1, ee + 6, e2);
75  CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(), &e2,
76  1, &e1);
77  }
78  EntityHandle f3, f4;
79  {
80  CHKERR m_field.get_moab().side_element(prism, 2, 3, f3);
81  CHKERR m_field.get_moab().side_element(prism, 2, 4, f4);
82  CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(), &f4,
83  1, &f3);
84  }
85  if (number_nodes > 3) {
86  EntityHandle meshset;
87  CHKERR m_field.get_moab().create_meshset(MESHSET_SET, meshset);
88  CHKERR m_field.get_moab().add_entities(meshset, &f4, 1);
89  for (int ee = 0; ee <= 2; ee++) {
90  EntityHandle e2;
91  CHKERR m_field.get_moab().side_element(prism, 1, ee + 6, e2);
92  CHKERR m_field.get_moab().add_entities(meshset, &e2, 1);
93  }
94  CHKERR m_field.get_moab().convert_entities(meshset, true, false, false);
95  CHKERR m_field.get_moab().delete_entities(&meshset, 1);
96  const EntityHandle *conn_f4;
97  int number_nodes_f4 = 0;
98  CHKERR m_field.get_moab().get_connectivity(f4, conn_f4, number_nodes_f4,
99  false);
100  if (number_nodes_f4 != number_nodes) {
101  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
102  "data inconsistency");
103  }
104  CHKERR m_field.get_moab().set_coords(&conn_f4[3], 3, &coords[9]);
105  }
106  }
107  MoFEMFunctionReturn(0);
108 }
109 
110 MoFEMErrorCode PrismsFromSurfaceInterface::seedPrismsEntities(
111  Range &prisms, const BitRefLevel &bit, int verb) {
112  Interface &m_field = cOre;
113  auto ref_ents_ptr = m_field.get_ref_ents();
114  MoFEMFunctionBegin;
115  MPI_Comm comm = m_field.get_comm();
116  RefEntity_multiIndex *refined_entities_ptr;
117  refined_entities_ptr =
118  const_cast<RefEntity_multiIndex *>(ref_ents_ptr);
119  if (!prisms.empty()) {
120  int dim = m_field.get_moab().dimension_from_handle(prisms[0]);
121  for (int dd = 0; dd <= dim; dd++) {
122  Range ents;
123  CHKERR m_field.get_moab().get_adjacencies(prisms, dd, true, ents,
124  moab::Interface::UNION);
125  Range::iterator eit = ents.begin();
126  for (; eit != ents.end(); eit++) {
127  std::pair<RefEntity_multiIndex::iterator, bool> p_ent =
128  refined_entities_ptr->insert(boost::shared_ptr<RefEntity>(
129  new RefEntity(m_field.get_basic_entity_data_ptr(), *eit)));
130  *(const_cast<RefEntity *>(p_ent.first->get())->getBitRefLevelPtr()) |=
131  bit;
132  if (verb >= VERY_VERBOSE) {
133  std::ostringstream ss;
134  ss << *(p_ent.first);
135  PetscSynchronizedPrintf(comm, "%s\n", ss.str().c_str());
136  }
137  }
138  }
139  }
140  MoFEMFunctionReturn(0);
141 }
142 
143 MoFEMErrorCode PrismsFromSurfaceInterface::createPrismsFromPrisms(
144  const Range &prisms, bool from_down, Range &out_prisms, int verb) {
145  MoFEMFunctionBegin;
146  Interface &m_field = cOre;
147  Range tris;
148  for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
149  EntityHandle face;
150  if (from_down) {
151  CHKERR m_field.get_moab().side_element(*pit, 2, 3, face);
152  } else {
153  CHKERR m_field.get_moab().side_element(*pit, 2, 4, face);
154  }
155  tris.insert(face);
156  }
157  CHKERR createPrisms(tris, NO_SWAP, out_prisms, verb);
158  MoFEMFunctionReturn(0);
159 }
160 
161 MoFEMErrorCode PrismsFromSurfaceInterface::setThickness(
162  const Range &prisms, const double director3[], const double director4[]) {
163  MoFEMFunctionBegin;
164  Interface &m_field = cOre;
165  Range nodes_f3, nodes_f4;
166  for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
167  for (int ff = 3; ff <= 4; ff++) {
168  EntityHandle face;
169  CHKERR m_field.get_moab().side_element(*pit, 2, ff, face);
170  const EntityHandle *conn;
171  int number_nodes = 0;
172  CHKERR m_field.get_moab().get_connectivity(face, conn, number_nodes,
173  false);
174  if (ff == 3) {
175  nodes_f3.insert(&conn[0], &conn[number_nodes]);
176  } else {
177  nodes_f4.insert(&conn[0], &conn[number_nodes]);
178  }
179  }
180  }
181  double coords[3];
182  for (Range::iterator nit = nodes_f3.begin(); nit != nodes_f3.end(); nit++) {
183  CHKERR m_field.get_moab().get_coords(&*nit, 1, coords);
184  cblas_daxpy(3, 1, director3, 1, coords, 1);
185  CHKERR m_field.get_moab().set_coords(&*nit, 1, coords);
186  }
187  for (Range::iterator nit = nodes_f4.begin(); nit != nodes_f4.end(); nit++) {
188  CHKERR m_field.get_moab().get_coords(&*nit, 1, coords);
189  cblas_daxpy(3, 1, director4, 1, coords, 1);
190  CHKERR m_field.get_moab().set_coords(&*nit, 1, coords);
191  }
192  MoFEMFunctionReturn(0);
193 }
194 
195 MoFEMErrorCode PrismsFromSurfaceInterface::setNormalThickness(
196  const Range &prisms, double thickness3, double thickness4) {
197  Interface &m_field = cOre;
198  MoFEMFunctionBegin;
199 
200  auto add_normal = [&](std::map<EntityHandle, std::array<double, 3>> &nodes,
201  EntityHandle face) {
202  MoFEMFunctionBeginHot;
203  const EntityHandle *conn;
204  int number_nodes;
205  CHKERR m_field.get_moab().get_connectivity(face, conn, number_nodes, false);
206  std::array<double, 9> coords;
207  CHKERR m_field.get_moab().get_coords(conn, number_nodes, coords.data());
208  std::array<double, 3> normal;
209  CHKERR Tools::getTriNormal(coords.data(), normal.data());
210  double a = sqrt(normal[0] * normal[0] + normal[1] * normal[1] +
211  normal[2] * normal[2]);
212  for (auto d : {0, 1, 2})
213  normal[d] /= a;
214  for (auto n : {0, 1, 2}) {
215  try {
216  for (auto d : {0, 1, 2})
217  nodes.at(conn[n])[d] += normal[d];
218  } catch (...) {
219  nodes.insert(
220  std::pair<EntityHandle, std::array<double, 3>>(conn[n], normal));
221  }
222  }
223  MoFEMFunctionReturnHot(0);
224  };
225 
226  auto apply_map = [&](auto &nodes, double t) {
227  MoFEMFunctionBeginHot;
228  for (auto &m : nodes) {
229  std::array<double, 3> coords;
230  CHKERR m_field.get_moab().get_coords(&m.first, 1, coords.data());
231  auto &normal = m.second;
232  double a = sqrt(normal[0] * normal[0] + normal[1] * normal[1] +
233  normal[2] * normal[2]);
234  for (auto d : {0, 1, 2})
235  coords[d] += (normal[d] / a) * t;
236  CHKERR m_field.get_moab().set_coords(&m.first, 1, coords.data());
237  }
238  MoFEMFunctionReturnHot(0);
239  };
240 
241  map<EntityHandle, std::array<double, 3>> nodes_f3, nodes_f4;
242  for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
243  for (int ff = 3; ff <= 4; ff++) {
244  EntityHandle face;
245  CHKERR m_field.get_moab().side_element(*pit, 2, ff, face);
246  if (ff == 3)
247  CHKERR add_normal(nodes_f3, face);
248  else
249  CHKERR add_normal(nodes_f4, face);
250  }
251  }
252 
253  CHKERR apply_map(nodes_f3, thickness3);
254  CHKERR apply_map(nodes_f4, thickness4);
255 
256  MoFEMFunctionReturn(0);
257 }
258 
259 MoFEMErrorCode
260 PrismsFromSurfaceInterface::updateMeshestByEdgeBlock(const Range &prisms) {
261  Interface &m_field = cOre;
262  MoFEMFunctionBegin;
263  Range prisms_edges;
264  CHKERR m_field.get_moab().get_adjacencies(prisms, 1, true, prisms_edges,
265  moab::Interface::UNION);
266  Range prisms_faces;
267  CHKERR m_field.get_moab().get_adjacencies(prisms, 2, true, prisms_faces,
268  moab::Interface::UNION);
269  for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, it)) {
270  Range edges;
271  CHKERR m_field.get_moab().get_entities_by_type(it->meshset, MBEDGE, edges,
272  true);
273  edges = intersect(edges, prisms_edges);
274  if (!edges.empty()) {
275  Range edges_faces;
276  CHKERR m_field.get_moab().get_adjacencies(edges, 2, false, edges_faces,
277  moab::Interface::UNION);
278  edges_faces = intersect(edges_faces, prisms_faces.subset_by_type(MBQUAD));
279  EntityHandle meshset = it->getMeshset();
280  CHKERR m_field.get_moab().add_entities(meshset, edges_faces);
281  }
282  }
283  MoFEMFunctionReturn(0);
284 }
285 
286 MoFEMErrorCode
287 PrismsFromSurfaceInterface::updateMeshestByTriBlock(const Range &prisms) {
288  Interface &m_field = cOre;
289  MoFEMFunctionBegin;
290  Range prisms_tris;
291  CHKERR m_field.get_moab().get_adjacencies(prisms, 2, true, prisms_tris,
292  moab::Interface::UNION);
293  prisms_tris = prisms_tris.subset_by_type(MBTRI);
294  for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, it)) {
295  Range tris;
296  CHKERR m_field.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
297  true);
298  tris = intersect(tris, prisms_tris);
299  if (!tris.empty()) {
300  Range tris_ents;
301  CHKERR m_field.get_moab().get_adjacencies(tris, 3, false, tris_ents,
302  moab::Interface::UNION);
303  tris_ents = intersect(tris_ents, prisms);
304  EntityHandle meshset = it->getMeshset();
305  CHKERR m_field.get_moab().add_entities(meshset, tris_ents);
306  }
307  }
308  MoFEMFunctionReturn(0);
309 }
310 
311 } // namespace MoFEM
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::PrismsFromSurfaceInterface
merge node from two bit levels
Definition: PrismsFromSurfaceInterface.hpp:15
DEPRECATED
#define DEPRECATED
Definition: definitions.h:17
EntityHandle
MoFEM::PrismsFromSurfaceInterface::SwapType
SwapType
List of types of node swapping performed on a created prism Node swapping is required to satisfy the ...
Definition: PrismsFromSurfaceInterface.hpp:37
MoFEM::PrismsFromSurfaceInterface::SWAP_TRI_NODE_ORDER
@ SWAP_TRI_NODE_ORDER
Definition: PrismsFromSurfaceInterface.hpp:39
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::PrismsFromSurfaceInterface::SWAP_TOP_AND_BOT_TRI
@ SWAP_TOP_AND_BOT_TRI
Definition: PrismsFromSurfaceInterface.hpp:40
MoFEM::RefEntityTmp< 0 >
Struct keeps handle to refined handle.
Definition: RefEntsMultiIndices.hpp:141
MoFEM::RefEntity_multiIndex
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_unique< tag< Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getEnt > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, ordered_non_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > > >, ordered_non_unique< tag< Composite_ParentEnt_And_EntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex
Definition: RefEntsMultiIndices.hpp:760
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::PrismsFromSurfaceInterface::updateMeshestByTriBlock
MoFEMErrorCode updateMeshestByTriBlock(const Range &prisms)
Add prism to bockset.
Definition: PrismsFromSurfaceInterface.cpp:287
MoFEM::PrismsFromSurfaceInterface::seedPrismsEntities
MoFEMErrorCode seedPrismsEntities(Range &prisms, const BitRefLevel &bit, int verb=-1)
Seed prism entities by bit level.
Definition: PrismsFromSurfaceInterface.cpp:110
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
MoFEM::PrismsFromSurfaceInterface::createPrismsFromPrisms
MoFEMErrorCode createPrismsFromPrisms(const Range &prisms, bool from_down, Range &out_prisms, int verb=-1)
Make prisms by extruding top or bottom prisms.
Definition: PrismsFromSurfaceInterface.cpp:143
MoFEM::CoreInterface::get_basic_entity_data_ptr
virtual boost::shared_ptr< BasicEntityData > & get_basic_entity_data_ptr()=0
Get pointer to basic entity data.
MoFEM::CoreTmp< 0 >::get_th_RefParentHandle
Tag get_th_RefParentHandle() const
Definition: Core.hpp:197
MoFEM::RefEntity
RefEntityTmp< 0 > RefEntity
Definition: RefEntsMultiIndices.hpp:566
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
MoFEM::PrismsFromSurfaceInterface::cOre
MoFEM::Core & cOre
Definition: PrismsFromSurfaceInterface.hpp:20
convert.n
n
Definition: convert.py:82
MoFEM::UnknownInterface
base class for all interface classes
Definition: UnknownInterface.hpp:34
Range
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
MoFEM::CoreInterface::get_ref_ents
virtual const RefEntity_multiIndex * get_ref_ents() const =0
Get the ref ents object.
MoFEM::PrismsFromSurfaceInterface::NO_SWAP
@ NO_SWAP
Definition: PrismsFromSurfaceInterface.hpp:38
MoFEM::PrismsFromSurfaceInterface::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition: PrismsFromSurfaceInterface.cpp:9
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
MoFEM::PrismsFromSurfaceInterface::createdVertices
std::map< EntityHandle, EntityHandle > createdVertices
Definition: PrismsFromSurfaceInterface.hpp:24
MoFEM::PrismsFromSurfaceInterface::updateMeshestByEdgeBlock
MoFEMErrorCode updateMeshestByEdgeBlock(const Range &prisms)
Add quads to bockset.
Definition: PrismsFromSurfaceInterface.cpp:260
_IT_CUBITMESHSETS_FOR_LOOP_
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
Definition: MeshsetsManager.hpp:34
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
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
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
MoFEM::PrismsFromSurfaceInterface::createPrisms
MoFEMErrorCode createPrisms(const Range &ents, const SwapType swap_type, Range &prisms, int verb=-1)
Make prisms from triangles.
Definition: PrismsFromSurfaceInterface.cpp:22
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEM::Tools::getTriNormal
static MoFEMErrorCode getTriNormal(const double *coords, double *normal, double *d_normal=nullptr)
Get the Tri Normal objectGet triangle normal.
Definition: Tools.cpp:353
VERY_VERBOSE
@ VERY_VERBOSE
Definition: definitions.h:223
MoFEM::PrismsFromSurfaceInterface::setThickness
MoFEMErrorCode setThickness(const Range &prisms, const double director3[], const double director4[])
Definition: PrismsFromSurfaceInterface.cpp:161
MoFEM::PrismsFromSurfaceInterface::setNormalThickness
MoFEMErrorCode setNormalThickness(const Range &prisms, double thickness3, double thickness4)
Definition: PrismsFromSurfaceInterface.cpp:195
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
Generated by Doxygen 1.8.17 and hosted at University of Glasgow