v0.9.1
PrismsFromSurfaceInterface.cpp
Go to the documentation of this file.
1 /** \file PrismsFromSurfaceInterface.cpp
2  * \brief Interface for creating prisms from surface elements
3  *
4  * MoFEM is free software: you can redistribute it and/or modify it under
5  * the terms of the GNU Lesser General Public License as published by the
6  * Free Software Foundation, either version 3 of the License, or (at your
7  * option) any later version.
8  *
9  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
12  * License for more details.
13  *
14  * You should have received a copy of the GNU Lesser General Public
15  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>
16  */
17 
18 namespace MoFEM {
19 
22  UnknownInterface **iface) const {
24  *iface = NULL;
25  if (uuid == IDD_MOFEMPrismsFromSurface) {
26  *iface = const_cast<PrismsFromSurfaceInterface *>(this);
28  }
29  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
31 }
32 
34  const Range &ents, Range &prisms, int verb) {
36  CHKERR createPrisms(ents, false, prisms, verb);
38 }
39 
41  const bool swap_nodes,
42  Range &prisms,
43  int verb) {
45  Interface &m_field = cOre;
46  Range tris = ents.subset_by_type(MBTRI);
47  for (Range::iterator tit = tris.begin(); tit != tris.end(); tit++) {
48  const EntityHandle *conn;
49  int number_nodes = 0;
50  CHKERR m_field.get_moab().get_connectivity(*tit, conn, number_nodes, false);
51  double coords[3 * number_nodes];
52  CHKERR m_field.get_moab().get_coords(conn, number_nodes, coords);
53  EntityHandle prism_nodes[6];
54  for (int nn = 0; nn < 3; nn++) {
55  prism_nodes[nn] = conn[nn];
56  if (createdVertices.find(conn[nn]) != createdVertices.end()) {
57  prism_nodes[3 + nn] = createdVertices[prism_nodes[nn]];
58  } else {
59  CHKERR m_field.get_moab().create_vertex(&coords[3 * nn],
60  prism_nodes[3 + nn]);
61  createdVertices[conn[nn]] = prism_nodes[3 + nn];
62  CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(),
63  &prism_nodes[3 + nn], 1,
64  &prism_nodes[nn]);
65  }
66  }
67  if (swap_nodes) {
68  std::swap(prism_nodes[1], prism_nodes[2]);
69  std::swap(prism_nodes[4], prism_nodes[5]);
70  }
71  EntityHandle prism;
72  CHKERR m_field.get_moab().create_element(MBPRISM, prism_nodes, 6, prism);
73  Range edges;
74  CHKERR m_field.get_moab().get_adjacencies(&prism, 1, 1, true, edges,
75  moab::Interface::UNION);
76  Range faces;
77  CHKERR m_field.get_moab().get_adjacencies(&prism, 1, 2, true, faces,
78  moab::Interface::UNION);
79  prisms.insert(prism);
80  for (int ee = 0; ee <= 2; ee++) {
81  EntityHandle e1;
82  CHKERR m_field.get_moab().side_element(prism, 1, ee, e1);
83  EntityHandle e2;
84  CHKERR m_field.get_moab().side_element(prism, 1, ee + 6, e2);
85  CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(), &e2,
86  1, &e1);
87  }
88  EntityHandle f3, f4;
89  {
90  CHKERR m_field.get_moab().side_element(prism, 2, 3, f3);
91  CHKERR m_field.get_moab().side_element(prism, 2, 4, f4);
92  CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(), &f4,
93  1, &f3);
94  }
95  if (number_nodes > 3) {
96  EntityHandle meshset;
97  CHKERR m_field.get_moab().create_meshset(
98  MESHSET_SET | MESHSET_TRACK_OWNER, meshset);
99  CHKERR m_field.get_moab().add_entities(meshset, &f4, 1);
100  for (int ee = 0; ee <= 2; ee++) {
101  EntityHandle e2;
102  CHKERR m_field.get_moab().side_element(prism, 1, ee + 6, e2);
103  CHKERR m_field.get_moab().add_entities(meshset, &e2, 1);
104  }
105  CHKERR m_field.get_moab().convert_entities(meshset, true, false, false);
106  CHKERR m_field.get_moab().delete_entities(&meshset, 1);
107  const EntityHandle *conn_f4;
108  int number_nodes_f4 = 0;
109  CHKERR m_field.get_moab().get_connectivity(f4, conn_f4, number_nodes_f4,
110  false);
111  if (number_nodes_f4 != number_nodes) {
112  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
113  "data inconsistency");
114  }
115  CHKERR m_field.get_moab().set_coords(&conn_f4[3], 3, &coords[9]);
116  }
117  }
119 }
120 
122  Range &prisms, const BitRefLevel &bit, int verb) {
124  Interface &m_field = cOre;
125  const RefEntity_multiIndex *const_refined_entities_ptr;
126  CHKERR m_field.get_ref_ents(&const_refined_entities_ptr);
127  MPI_Comm comm = m_field.get_comm();
128  RefEntity_multiIndex *refined_entities_ptr;
129  refined_entities_ptr =
130  const_cast<RefEntity_multiIndex *>(const_refined_entities_ptr);
131  if (!prisms.empty()) {
132  int dim = m_field.get_moab().dimension_from_handle(prisms[0]);
133  for (int dd = 0; dd <= dim; dd++) {
134  Range ents;
135  CHKERR m_field.get_moab().get_adjacencies(prisms, dd, true, ents,
136  moab::Interface::UNION);
137  Range::iterator eit = ents.begin();
138  for (; eit != ents.end(); eit++) {
139  std::pair<RefEntity_multiIndex::iterator, bool> p_ent =
140  refined_entities_ptr->insert(boost::shared_ptr<RefEntity>(
141  new RefEntity(m_field.get_basic_entity_data_ptr(), *eit)));
142  *(const_cast<RefEntity *>(p_ent.first->get())->getBitRefLevelPtr()) |=
143  bit;
144  if (verb >= VERY_VERBOSE) {
145  std::ostringstream ss;
146  ss << *(p_ent.first);
147  PetscSynchronizedPrintf(comm, "%s\n", ss.str().c_str());
148  }
149  }
150  }
151  }
153 }
154 
156  const Range &prisms, bool from_down, Range &out_prisms, int verb) {
158  Interface &m_field = cOre;
159  Range tris;
160  for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
161  EntityHandle face;
162  if (from_down) {
163  CHKERR m_field.get_moab().side_element(*pit, 2, 3, face);
164  } else {
165  CHKERR m_field.get_moab().side_element(*pit, 2, 4, face);
166  }
167  tris.insert(face);
168  }
169  CHKERR createPrisms(tris, out_prisms, verb);
171 }
172 
174  const Range &prisms, const double director3[], const double director4[]) {
176  Interface &m_field = cOre;
177  Range nodes_f3, nodes_f4;
178  for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
179  for (int ff = 3; ff <= 4; ff++) {
180  EntityHandle face;
181  CHKERR m_field.get_moab().side_element(*pit, 2, ff, face);
182  const EntityHandle *conn;
183  int number_nodes = 0;
184  CHKERR m_field.get_moab().get_connectivity(face, conn, number_nodes,
185  false);
186  if (ff == 3) {
187  nodes_f3.insert(&conn[0], &conn[number_nodes]);
188  } else {
189  nodes_f4.insert(&conn[0], &conn[number_nodes]);
190  }
191  }
192  }
193  double coords[3];
194  for (Range::iterator nit = nodes_f3.begin(); nit != nodes_f3.end(); nit++) {
195  CHKERR m_field.get_moab().get_coords(&*nit, 1, coords);
196  cblas_daxpy(3, 1, director3, 1, coords, 1);
197  CHKERR m_field.get_moab().set_coords(&*nit, 1, coords);
198  }
199  for (Range::iterator nit = nodes_f4.begin(); nit != nodes_f4.end(); nit++) {
200  CHKERR m_field.get_moab().get_coords(&*nit, 1, coords);
201  cblas_daxpy(3, 1, director4, 1, coords, 1);
202  CHKERR m_field.get_moab().set_coords(&*nit, 1, coords);
203  }
205 }
206 
208  const Range &prisms, double thickness3, double thickness4) {
209  Interface &m_field = cOre;
211 
212  auto add_normal = [&](std::map<EntityHandle, std::array<double, 3>> &nodes,
213  EntityHandle face) {
215  const EntityHandle *conn;
216  int number_nodes;
217  CHKERR m_field.get_moab().get_connectivity(face, conn, number_nodes, false);
218  std::array<double, 9> coords;
219  CHKERR m_field.get_moab().get_coords(conn, number_nodes, coords.data());
220  std::array<double, 3> normal;
221  CHKERR Tools::getTriNormal(coords.data(), normal.data());
222  double a = sqrt(normal[0] * normal[0] + normal[1] * normal[1] +
223  normal[2] * normal[2]);
224  for (auto d : {0, 1, 2})
225  normal[d] /= a;
226  for (auto n : {0, 1, 2}) {
227  try {
228  for (auto d : {0, 1, 2})
229  nodes.at(conn[n])[d] += normal[d];
230  } catch (...) {
231  nodes.insert(
232  std::pair<EntityHandle, std::array<double, 3>>(conn[n], normal));
233  }
234  }
236  };
237 
238  auto apply_map = [&](auto &nodes, double t) {
240  for (auto &m : nodes) {
241  std::array<double, 3> coords;
242  CHKERR m_field.get_moab().get_coords(&m.first, 1, coords.data());
243  auto &normal = m.second;
244  double a = sqrt(normal[0] * normal[0] + normal[1] * normal[1] +
245  normal[2] * normal[2]);
246  for (auto d : {0, 1, 2})
247  coords[d] += (normal[d] / a) * t;
248  CHKERR m_field.get_moab().set_coords(&m.first, 1, coords.data());
249  }
251  };
252 
253  map<EntityHandle, std::array<double, 3>> nodes_f3, nodes_f4;
254  for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
255  for (int ff = 3; ff <= 4; ff++) {
256  EntityHandle face;
257  CHKERR m_field.get_moab().side_element(*pit, 2, ff, face);
258  if (ff == 3)
259  CHKERR add_normal(nodes_f3, face);
260  else
261  CHKERR add_normal(nodes_f4, face);
262  }
263  }
264 
265  CHKERR apply_map(nodes_f3, thickness3);
266  CHKERR apply_map(nodes_f4, thickness4);
267 
269 }
270 
273  Interface &m_field = cOre;
275  Range prisms_edges;
276  CHKERR m_field.get_moab().get_adjacencies(prisms, 1, true, prisms_edges,
277  moab::Interface::UNION);
278  Range prisms_faces;
279  CHKERR m_field.get_moab().get_adjacencies(prisms, 2, true, prisms_faces,
280  moab::Interface::UNION);
281  for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, it)) {
282  Range edges;
283  CHKERR m_field.get_moab().get_entities_by_type(it->meshset, MBEDGE, edges,
284  true);
285  edges = intersect(edges, prisms_edges);
286  if (!edges.empty()) {
287  Range edges_faces;
288  CHKERR m_field.get_moab().get_adjacencies(edges, 2, false, edges_faces,
289  moab::Interface::UNION);
290  edges_faces = intersect(edges_faces, prisms_faces.subset_by_type(MBQUAD));
291  EntityHandle meshset = it->getMeshset();
292  CHKERR m_field.get_moab().add_entities(meshset, edges_faces);
293  }
294  }
296 }
297 
300  Interface &m_field = cOre;
302  Range prisms_tris;
303  CHKERR m_field.get_moab().get_adjacencies(prisms, 2, true, prisms_tris,
304  moab::Interface::UNION);
305  prisms_tris = prisms_tris.subset_by_type(MBTRI);
306  for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, it)) {
307  Range tris;
308  CHKERR m_field.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
309  true);
310  tris = intersect(tris, prisms_tris);
311  if (!tris.empty()) {
312  Range tris_ents;
313  CHKERR m_field.get_moab().get_adjacencies(tris, 3, false, tris_ents,
314  moab::Interface::UNION);
315  tris_ents = intersect(tris_ents, prisms);
316  EntityHandle meshset = it->getMeshset();
317  CHKERR m_field.get_moab().add_entities(meshset, tris_ents);
318  }
319  }
321 }
322 
323 } // namespace MoFEM
virtual MoFEMErrorCode get_ref_ents(const RefEntity_multiIndex **refined_ents_ptr) const =0
Get ref entities multi-index from database.
Deprecated interface functions.
MoFEM interface unique ID.
Tag get_th_RefParentHandle() const
Definition: Core.hpp:150
virtual moab::Interface & get_moab()=0
virtual boost::shared_ptr< BasicEntityData > & get_basic_entity_data_ptr()=0
Get pointer to basic entity data.
MoFEMErrorCode setThickness(const Range &prisms, const double director3[], const double director4[])
void cblas_daxpy(const int N, const double alpha, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_daxpy.c:11
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:507
MoFEMErrorCode seedPrismsEntities(Range &prisms, const BitRefLevel &bit, int verb=-1)
Seed prism entities by bit level.
base class for all interface classes
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
FTensor::Tensor1< double, 2 > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, 2 > &t_disp)
Definition: ContactOps.hpp:128
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
MoFEMErrorCode updateMeshestByTriBlock(const Range &prisms)
Add prism to bockset.
static const MOFEMuuid IDD_MOFEMPrismsFromSurface
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
Struct keeps handle to refined handle.
MoFEMErrorCode createPrisms(const Range &ents, const bool swap_nodes, Range &prisms, int verb=-1)
Make prisms from triangles.
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MoFEMErrorCode setNormalThickness(const Range &prisms, double thickness3, double thickness4)
const int dim
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
MoFEMErrorCode updateMeshestByEdgeBlock(const Range &prisms)
Add quads to bockset.
MoFEMErrorCode query_interface(const MOFEMuuid &uuid, UnknownInterface **iface) const
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
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
std::map< EntityHandle, EntityHandle > createdVertices
static MoFEMErrorCode getTriNormal(const double *coords, double *normal)
Get the Tri Normal objectGet triangle normal.
Definition: Tools.cpp:258
MoFEMErrorCode createPrismsFromPrisms(const Range &prisms, bool from_down, Range &out_prisms, int verb=-1)
Make prisms by extruding top or bottom prisms.
#define CHKERR
Inline error check.
Definition: definitions.h:602
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:50
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_unique< tag< Ent_mi_tag >, member< RefEntity::BasicEntity, EntityHandle, &RefEntity::ent > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType > >, ordered_non_unique< tag< ParentEntType_mi_tag >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > >, ordered_non_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, 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::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex
#define DEPRECATED
Definition: definitions.h:29
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
virtual MPI_Comm & get_comm() const =0