v0.9.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  * 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  Range &prisms,
35  int verb) {
37  Interface &m_field = cOre;
38  Range tris = ents.subset_by_type(MBTRI);
39  for (Range::iterator tit = tris.begin(); tit != tris.end(); tit++) {
40  const EntityHandle *conn;
41  int number_nodes = 0;
42  CHKERR m_field.get_moab().get_connectivity(*tit, conn, number_nodes, false);
43  double coords[3 * number_nodes];
44  CHKERR m_field.get_moab().get_coords(conn, number_nodes, coords);
45  EntityHandle prism_nodes[6];
46  for (int nn = 0; nn < 3; nn++) {
47  prism_nodes[nn] = conn[nn];
48  if (createdVertices.find(conn[nn]) != createdVertices.end()) {
49  prism_nodes[3 + nn] = createdVertices[prism_nodes[nn]];
50  } else {
51  CHKERR m_field.get_moab().create_vertex(&coords[3 * nn],
52  prism_nodes[3 + nn]);
53  createdVertices[conn[nn]] = prism_nodes[3 + nn];
54  CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(),
55  &prism_nodes[3 + nn], 1,
56  &prism_nodes[nn]);
57  }
58  }
59  EntityHandle prism;
60  CHKERR m_field.get_moab().create_element(MBPRISM, prism_nodes, 6, prism);
61  Range edges;
62  CHKERR m_field.get_moab().get_adjacencies(&prism, 1, 1, true, edges,
63  moab::Interface::UNION);
64  Range faces;
65  CHKERR m_field.get_moab().get_adjacencies(&prism, 1, 2, true, faces,
66  moab::Interface::UNION);
67  prisms.insert(prism);
68  for (int ee = 0; ee <= 2; ee++) {
69  EntityHandle e1;
70  CHKERR m_field.get_moab().side_element(prism, 1, ee, e1);
71  EntityHandle e2;
72  CHKERR m_field.get_moab().side_element(prism, 1, ee + 6, e2);
73  CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(), &e2,
74  1, &e1);
75  }
76  EntityHandle f3, f4;
77  {
78  CHKERR m_field.get_moab().side_element(prism, 2, 3, f3);
79  CHKERR m_field.get_moab().side_element(prism, 2, 4, f4);
80  CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(), &f4,
81  1, &f3);
82  }
83  if (number_nodes > 3) {
84  EntityHandle meshset;
85  CHKERR m_field.get_moab().create_meshset(
86  MESHSET_SET | MESHSET_TRACK_OWNER, meshset);
87  CHKERR m_field.get_moab().add_entities(meshset, &f4, 1);
88  for (int ee = 0; ee <= 2; ee++) {
89  EntityHandle e2;
90  CHKERR m_field.get_moab().side_element(prism, 1, ee + 6, e2);
91  CHKERR m_field.get_moab().add_entities(meshset, &e2, 1);
92  }
93  CHKERR m_field.get_moab().convert_entities(meshset, true, false, false);
94  CHKERR m_field.get_moab().delete_entities(&meshset, 1);
95  const EntityHandle *conn_f4;
96  int number_nodes_f4 = 0;
97  CHKERR m_field.get_moab().get_connectivity(f4, conn_f4, number_nodes_f4,
98  false);
99  if (number_nodes_f4 != number_nodes) {
100  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
101  "data inconsistency");
102  }
103  CHKERR m_field.get_moab().set_coords(&conn_f4[3], 3, &coords[9]);
104  }
105  }
107 }
108 
110  Range &prisms, const BitRefLevel &bit, int verb) {
112  Interface &m_field = cOre;
113  const RefEntity_multiIndex *const_refined_entities_ptr;
114  CHKERR m_field.get_ref_ents(&const_refined_entities_ptr);
115  MPI_Comm comm = m_field.get_comm();
116  RefEntity_multiIndex *refined_entities_ptr;
117  refined_entities_ptr =
118  const_cast<RefEntity_multiIndex *>(const_refined_entities_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  }
141 }
142 
144  const Range &prisms, bool from_down, Range &out_prisms, int verb) {
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, out_prisms, verb);
159 }
160 
162  const Range &prisms, const double director3[], const double director4[]) {
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  }
193 }
194 
196  const Range &prisms, double thickness3, double thickness4) {
197  Interface &m_field = cOre;
199 
200  auto add_normal = [&](std::map<EntityHandle, std::array<double, 3>> &nodes,
201  EntityHandle face) {
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  }
224  };
225 
226  auto apply_map = [&](auto &nodes, double t) {
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  }
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 
257 }
258 
261  Interface &m_field = cOre;
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  }
284 }
285 
288  Interface &m_field = cOre;
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  }
309 }
310 
311 } // namespace MoFEM
virtual MoFEMErrorCode get_ref_ents(const RefEntity_multiIndex **refined_ents_ptr) const =0
Get ref entities multi-index from database.
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:501
MoFEMErrorCode seedPrismsEntities(Range &prisms, const BitRefLevel &bit, int verb=-1)
Seed prism entities by bit level.
base class for all interface classes
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
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:508
Struct keeps handle to refined handle.
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MoFEMErrorCode setNormalThickness(const Range &prisms, double thickness3, double thickness4)
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
MoFEMErrorCode createPrisms(const Range &ents, Range &prisms, int verb=-1)
Make prisms from triangles.
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:596
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
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 _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:407
virtual MPI_Comm & get_comm() const =0