v0.13.1
Loading...
Searching...
No Matches
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
7namespace MoFEM {
8
10 boost::typeindex::type_index type_index, UnknownInterface **iface) const {
11 *iface = const_cast<PrismsFromSurfaceInterface *>(this);
12 return 0;
13}
14
16 const Range &ents, Range &prisms, int verb) {
18 CHKERR createPrisms(ents, NO_SWAP, prisms, verb);
20}
21
23 const Range &ents, const SwapType swap_type, Range &prisms, int verb) {
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) {
49 std::swap(prism_nodes[1], prism_nodes[2]);
50 std::swap(prism_nodes[4], prism_nodes[5]);
51 break;
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 }
108}
109
111 Range &prisms, const BitRefLevel &bit, int verb) {
112 Interface &m_field = cOre;
113 auto ref_ents_ptr = m_field.get_ref_ents();
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 }
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
constexpr double a
@ VERY_VERBOSE
Definition: definitions.h:210
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#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 DEPRECATED
Definition: definitions.h:17
#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
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
const int dim
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
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
virtual const RefEntity_multiIndex * get_ref_ents() const =0
Get the ref ents object.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
auto bit
set bit
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
RefEntityTmp< 0 > RefEntity
constexpr double t
plate stiffness
Definition: plate.cpp:59
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual boost::shared_ptr< BasicEntityData > & get_basic_entity_data_ptr()=0
Get pointer to basic entity data.
Tag get_th_RefParentHandle() const
Definition: Core.hpp:197
Deprecated interface functions.
std::map< EntityHandle, EntityHandle > createdVertices
MoFEMErrorCode setNormalThickness(const Range &prisms, double thickness3, double thickness4)
MoFEMErrorCode createPrisms(const Range &ents, const SwapType swap_type, Range &prisms, int verb=-1)
Make prisms from triangles.
MoFEMErrorCode updateMeshestByEdgeBlock(const Range &prisms)
Add quads to bockset.
MoFEMErrorCode setThickness(const Range &prisms, const double director3[], const double director4[])
SwapType
List of types of node swapping performed on a created prism Node swapping is required to satisfy the ...
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
MoFEMErrorCode seedPrismsEntities(Range &prisms, const BitRefLevel &bit, int verb=-1)
Seed prism entities by bit level.
MoFEMErrorCode createPrismsFromPrisms(const Range &prisms, bool from_down, Range &out_prisms, int verb=-1)
Make prisms by extruding top or bottom prisms.
MoFEMErrorCode updateMeshestByTriBlock(const Range &prisms)
Add prism to bockset.
Struct keeps handle to refined handle.
static MoFEMErrorCode getTriNormal(const double *coords, double *normal)
Get the Tri Normal objectGet triangle normal.
Definition: Tools.cpp:352
base class for all interface classes