v0.13.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
18namespace MoFEM {
19
21 boost::typeindex::type_index type_index, UnknownInterface **iface) const {
22 *iface = const_cast<PrismsFromSurfaceInterface *>(this);
23 return 0;
24}
25
27 const Range &ents, Range &prisms, int verb) {
29 CHKERR createPrisms(ents, NO_SWAP, prisms, verb);
31}
32
34 const Range &ents, const SwapType swap_type, Range &prisms, int verb) {
36 Interface &m_field = cOre;
37 Range tris = ents.subset_by_type(MBTRI);
38 for (Range::iterator tit = tris.begin(); tit != tris.end(); tit++) {
39 const EntityHandle *conn;
40 int number_nodes = 0;
41 CHKERR m_field.get_moab().get_connectivity(*tit, conn, number_nodes, false);
42 double coords[3 * number_nodes];
43 CHKERR m_field.get_moab().get_coords(conn, number_nodes, coords);
44 EntityHandle prism_nodes[6];
45 for (int nn = 0; nn < 3; nn++) {
46 prism_nodes[nn] = conn[nn];
47 if (createdVertices.find(conn[nn]) != createdVertices.end()) {
48 prism_nodes[3 + nn] = createdVertices[prism_nodes[nn]];
49 } else {
50 CHKERR m_field.get_moab().create_vertex(&coords[3 * nn],
51 prism_nodes[3 + nn]);
52 createdVertices[conn[nn]] = prism_nodes[3 + nn];
53 CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(),
54 &prism_nodes[3 + nn], 1,
55 &prism_nodes[nn]);
56 }
57 }
58 switch (swap_type) {
60 std::swap(prism_nodes[1], prism_nodes[2]);
61 std::swap(prism_nodes[4], prism_nodes[5]);
62 break;
64 std::swap(prism_nodes[0], prism_nodes[3]);
65 std::swap(prism_nodes[1], prism_nodes[4]);
66 std::swap(prism_nodes[2], prism_nodes[5]);
67 break;
68 case NO_SWAP:
69 default:
70 break;
71 }
72 EntityHandle prism;
73 CHKERR m_field.get_moab().create_element(MBPRISM, prism_nodes, 6, prism);
74 Range edges;
75 CHKERR m_field.get_moab().get_adjacencies(&prism, 1, 1, true, edges,
76 moab::Interface::UNION);
77 Range faces;
78 CHKERR m_field.get_moab().get_adjacencies(&prism, 1, 2, true, faces,
79 moab::Interface::UNION);
80 prisms.insert(prism);
81 for (int ee = 0; ee <= 2; ee++) {
82 EntityHandle e1;
83 CHKERR m_field.get_moab().side_element(prism, 1, ee, e1);
84 EntityHandle e2;
85 CHKERR m_field.get_moab().side_element(prism, 1, ee + 6, e2);
86 CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(), &e2,
87 1, &e1);
88 }
89 EntityHandle f3, f4;
90 {
91 CHKERR m_field.get_moab().side_element(prism, 2, 3, f3);
92 CHKERR m_field.get_moab().side_element(prism, 2, 4, f4);
93 CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(), &f4,
94 1, &f3);
95 }
96 if (number_nodes > 3) {
97 EntityHandle meshset;
98 CHKERR m_field.get_moab().create_meshset(MESHSET_SET, 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) {
123 Interface &m_field = cOre;
124 auto ref_ents_ptr = m_field.get_ref_ents();
126 MPI_Comm comm = m_field.get_comm();
127 RefEntity_multiIndex *refined_entities_ptr;
128 refined_entities_ptr =
129 const_cast<RefEntity_multiIndex *>(ref_ents_ptr);
130 if (!prisms.empty()) {
131 int dim = m_field.get_moab().dimension_from_handle(prisms[0]);
132 for (int dd = 0; dd <= dim; dd++) {
133 Range ents;
134 CHKERR m_field.get_moab().get_adjacencies(prisms, dd, true, ents,
135 moab::Interface::UNION);
136 Range::iterator eit = ents.begin();
137 for (; eit != ents.end(); eit++) {
138 std::pair<RefEntity_multiIndex::iterator, bool> p_ent =
139 refined_entities_ptr->insert(boost::shared_ptr<RefEntity>(
140 new RefEntity(m_field.get_basic_entity_data_ptr(), *eit)));
141 *(const_cast<RefEntity *>(p_ent.first->get())->getBitRefLevelPtr()) |=
142 bit;
143 if (verb >= VERY_VERBOSE) {
144 std::ostringstream ss;
145 ss << *(p_ent.first);
146 PetscSynchronizedPrintf(comm, "%s\n", ss.str().c_str());
147 }
148 }
149 }
150 }
152}
153
155 const Range &prisms, bool from_down, Range &out_prisms, int verb) {
157 Interface &m_field = cOre;
158 Range tris;
159 for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
160 EntityHandle face;
161 if (from_down) {
162 CHKERR m_field.get_moab().side_element(*pit, 2, 3, face);
163 } else {
164 CHKERR m_field.get_moab().side_element(*pit, 2, 4, face);
165 }
166 tris.insert(face);
167 }
168 CHKERR createPrisms(tris, out_prisms, verb);
170}
171
173 const Range &prisms, const double director3[], const double director4[]) {
175 Interface &m_field = cOre;
176 Range nodes_f3, nodes_f4;
177 for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
178 for (int ff = 3; ff <= 4; ff++) {
179 EntityHandle face;
180 CHKERR m_field.get_moab().side_element(*pit, 2, ff, face);
181 const EntityHandle *conn;
182 int number_nodes = 0;
183 CHKERR m_field.get_moab().get_connectivity(face, conn, number_nodes,
184 false);
185 if (ff == 3) {
186 nodes_f3.insert(&conn[0], &conn[number_nodes]);
187 } else {
188 nodes_f4.insert(&conn[0], &conn[number_nodes]);
189 }
190 }
191 }
192 double coords[3];
193 for (Range::iterator nit = nodes_f3.begin(); nit != nodes_f3.end(); nit++) {
194 CHKERR m_field.get_moab().get_coords(&*nit, 1, coords);
195 cblas_daxpy(3, 1, director3, 1, coords, 1);
196 CHKERR m_field.get_moab().set_coords(&*nit, 1, coords);
197 }
198 for (Range::iterator nit = nodes_f4.begin(); nit != nodes_f4.end(); nit++) {
199 CHKERR m_field.get_moab().get_coords(&*nit, 1, coords);
200 cblas_daxpy(3, 1, director4, 1, coords, 1);
201 CHKERR m_field.get_moab().set_coords(&*nit, 1, coords);
202 }
204}
205
207 const Range &prisms, double thickness3, double thickness4) {
208 Interface &m_field = cOre;
210
211 auto add_normal = [&](std::map<EntityHandle, std::array<double, 3>> &nodes,
212 EntityHandle face) {
214 const EntityHandle *conn;
215 int number_nodes;
216 CHKERR m_field.get_moab().get_connectivity(face, conn, number_nodes, false);
217 std::array<double, 9> coords;
218 CHKERR m_field.get_moab().get_coords(conn, number_nodes, coords.data());
219 std::array<double, 3> normal;
220 CHKERR Tools::getTriNormal(coords.data(), normal.data());
221 double a = sqrt(normal[0] * normal[0] + normal[1] * normal[1] +
222 normal[2] * normal[2]);
223 for (auto d : {0, 1, 2})
224 normal[d] /= a;
225 for (auto n : {0, 1, 2}) {
226 try {
227 for (auto d : {0, 1, 2})
228 nodes.at(conn[n])[d] += normal[d];
229 } catch (...) {
230 nodes.insert(
231 std::pair<EntityHandle, std::array<double, 3>>(conn[n], normal));
232 }
233 }
235 };
236
237 auto apply_map = [&](auto &nodes, double t) {
239 for (auto &m : nodes) {
240 std::array<double, 3> coords;
241 CHKERR m_field.get_moab().get_coords(&m.first, 1, coords.data());
242 auto &normal = m.second;
243 double a = sqrt(normal[0] * normal[0] + normal[1] * normal[1] +
244 normal[2] * normal[2]);
245 for (auto d : {0, 1, 2})
246 coords[d] += (normal[d] / a) * t;
247 CHKERR m_field.get_moab().set_coords(&m.first, 1, coords.data());
248 }
250 };
251
252 map<EntityHandle, std::array<double, 3>> nodes_f3, nodes_f4;
253 for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
254 for (int ff = 3; ff <= 4; ff++) {
255 EntityHandle face;
256 CHKERR m_field.get_moab().side_element(*pit, 2, ff, face);
257 if (ff == 3)
258 CHKERR add_normal(nodes_f3, face);
259 else
260 CHKERR add_normal(nodes_f4, face);
261 }
262 }
263
264 CHKERR apply_map(nodes_f3, thickness3);
265 CHKERR apply_map(nodes_f4, thickness4);
266
268}
269
272 Interface &m_field = cOre;
274 Range prisms_edges;
275 CHKERR m_field.get_moab().get_adjacencies(prisms, 1, true, prisms_edges,
276 moab::Interface::UNION);
277 Range prisms_faces;
278 CHKERR m_field.get_moab().get_adjacencies(prisms, 2, true, prisms_faces,
279 moab::Interface::UNION);
280 for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, it)) {
281 Range edges;
282 CHKERR m_field.get_moab().get_entities_by_type(it->meshset, MBEDGE, edges,
283 true);
284 edges = intersect(edges, prisms_edges);
285 if (!edges.empty()) {
286 Range edges_faces;
287 CHKERR m_field.get_moab().get_adjacencies(edges, 2, false, edges_faces,
288 moab::Interface::UNION);
289 edges_faces = intersect(edges_faces, prisms_faces.subset_by_type(MBQUAD));
290 EntityHandle meshset = it->getMeshset();
291 CHKERR m_field.get_moab().add_entities(meshset, edges_faces);
292 }
293 }
295}
296
299 Interface &m_field = cOre;
301 Range prisms_tris;
302 CHKERR m_field.get_moab().get_adjacencies(prisms, 2, true, prisms_tris,
303 moab::Interface::UNION);
304 prisms_tris = prisms_tris.subset_by_type(MBTRI);
305 for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, it)) {
306 Range tris;
307 CHKERR m_field.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
308 true);
309 tris = intersect(tris, prisms_tris);
310 if (!tris.empty()) {
311 Range tris_ents;
312 CHKERR m_field.get_moab().get_adjacencies(tris, 3, false, tris_ents,
313 moab::Interface::UNION);
314 tris_ents = intersect(tris_ents, prisms);
315 EntityHandle meshset = it->getMeshset();
316 CHKERR m_field.get_moab().add_entities(meshset, tris_ents);
317 }
318 }
320}
321
322} // namespace MoFEM
static Index< 'n', 3 > n
constexpr double a
@ VERY_VERBOSE
Definition: definitions.h:223
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#define DEPRECATED
Definition: definitions.h:30
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
const int dim
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
FTensor::Tensor1< double, SPACE_DIM > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, SPACE_DIM > &t_disp)
Definition: ContactOps.hpp:77
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
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:67
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
RefEntityTmp< 0 > RefEntity
constexpr double t
plate stiffness
Definition: plate.cpp:76
FTensor::Index< 'm', 3 > m
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:207
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:363
base class for all interface classes