v0.9.0
Classes | Typedefs | Functions | Variables
prism_elements_from_surface.cpp File Reference
#include <MoFEM.hpp>

Go to the source code of this file.

Classes

struct  CoordsAndHandle
 
struct  PrismOp
 
struct  PrismFE
 
struct  Op< OP >
 
struct  TriFE
 
struct  QuadFE
 
struct  EdgeFE
 

Typedefs

typedef multi_index_container< CoordsAndHandle, indexed_by< hashed_unique< composite_key< CoordsAndHandle, member< CoordsAndHandle, int, &CoordsAndHandle::x >, member< CoordsAndHandle, int, &CoordsAndHandle::y >, member< CoordsAndHandle, int, &CoordsAndHandle::z > > > > > MapCoords
 

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 
static constexpr int precision_exponent = 5
 
static constexpr int number_of_prisms_layers = 18
 
static constexpr double delta
 
static constexpr std::array< double, 3 > d3 = {0, 0, 0}
 
static constexpr std::array< double, 3 > d4 = {0, 0, delta}
 

Typedef Documentation

◆ MapCoords

typedef multi_index_container< CoordsAndHandle, indexed_by< hashed_unique<composite_key< CoordsAndHandle, member<CoordsAndHandle, int, &CoordsAndHandle::x>, member<CoordsAndHandle, int, &CoordsAndHandle::y>, member<CoordsAndHandle, int, &CoordsAndHandle::z> > > > > MapCoords

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)
Examples
prism_elements_from_surface.cpp.

Definition at line 143 of file prism_elements_from_surface.cpp.

143  {
144 
145  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
146 
147  try {
148 
149  moab::Core mb_instance;
150  moab::Interface &moab = mb_instance;
151  int rank;
152  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
153 
154  // Read parameters from line command
155  PetscBool flg = PETSC_TRUE;
156  char mesh_file_name[255];
157  CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
158  255, &flg);
159  if (flg != PETSC_TRUE)
160  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
161  "error -my_file (MESH FILE NEEDED)");
162 
163  // Read mesh to MOAB
164  const char *option;
165  option = ""; //"PARALLEL=BCAST;";//;DEBUG_IO";
166  CHKERR moab.load_file(mesh_file_name, 0, option);
167  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
168  if (pcomm == NULL)
169  pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
170 
171  Range tris;
172  CHKERR moab.get_entities_by_type(0, MBTRI, tris, false);
173  Range tris_verts;
174  CHKERR moab.get_connectivity(tris, tris_verts);
175  MatrixDouble tri_coords(tris_verts.size(), 3);
176  CHKERR moab.get_coords(tris_verts, &tri_coords(0, 0));
177 
178  MoFEM::Core core(moab);
179  MoFEM::Interface &m_field = core;
180 
181  std::array<Range, 3> edge_block;
182  for (auto b : {1, 2, 3})
184  b, BLOCKSET, 1, edge_block[b - 1]);
185 
186  PrismsFromSurfaceInterface *prisms_from_surface_interface;
187  CHKERR m_field.getInterface(prisms_from_surface_interface);
188 
189  Range prisms;
190  CHKERR prisms_from_surface_interface->createPrisms(tris, prisms);
191  prisms_from_surface_interface->setThickness(prisms, d3.data(), d4.data());
192  Range add_prims_layer;
193  Range extrude_prisms = prisms;
194 
195  for (int ll = 1; ll != number_of_prisms_layers; ++ll) {
196  prisms_from_surface_interface->createdVertices.clear();
197  CHKERR prisms_from_surface_interface->createPrismsFromPrisms(
198  extrude_prisms, false, add_prims_layer);
199  prisms_from_surface_interface->setThickness(add_prims_layer, d3.data(),
200  d4.data());
201  extrude_prisms = add_prims_layer;
202  prisms.merge(add_prims_layer);
203  add_prims_layer.clear();
204  }
205 
206  EntityHandle meshset;
207  CHKERR moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, meshset);
208  CHKERR moab.add_entities(meshset, prisms);
209 
210  MapCoords map_coords;
211  Range verts;
212  CHKERR moab.get_connectivity(prisms, verts);
213  MatrixDouble coords(verts.size(), 3);
214  CHKERR moab.get_coords(verts, &coords(0, 0));
215 
216  for (size_t v = 0; v != verts.size(); ++v)
217  map_coords.insert(CoordsAndHandle(&coords(v, 0), verts[v]));
218 
219  EntityHandle one_prism_meshset;
220  CHKERR moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER,
221  one_prism_meshset);
222 
223  std::array<double, 18> one_prism_coords = {0, 0, 0, 1, 0, 0, 0, 1, 0,
224  0, 0, 1, 1, 0, 1, 0, 1, 1};
225  std::array<EntityHandle, 6> one_prism_nodes;
226  for (int n = 0; n != 6; ++n)
227  CHKERR moab.create_vertex(&one_prism_coords[3 * n], one_prism_nodes[n]);
228  EntityHandle one_prism;
229  CHKERR m_field.get_moab().create_element(MBPRISM, one_prism_nodes.data(), 6,
230  one_prism);
231  Range one_prism_range;
232  one_prism_range.insert(one_prism);
233  CHKERR moab.add_entities(one_prism_meshset, one_prism_range);
234  Range one_prism_verts;
235  CHKERR moab.get_connectivity(one_prism_range, one_prism_verts);
236  CHKERR moab.add_entities(one_prism_meshset, one_prism_verts);
237  Range one_prism_adj_ents;
238  for (int d = 1; d != 3; ++d)
239  CHKERR moab.get_adjacencies(one_prism_range, d, true, one_prism_adj_ents,
240  moab::Interface::UNION);
241  CHKERR moab.add_entities(one_prism_meshset, one_prism_adj_ents);
242 
243  BitRefLevel bit_level0;
244  bit_level0.set(0);
245  CHKERR m_field.getInterface<BitRefManager>()->setEntitiesBitRefLevel(
246  one_prism_range, bit_level0);
247  CHKERR prisms_from_surface_interface->seedPrismsEntities(one_prism_range,
248  bit_level0);
249 
250  // Fields
251  CHKERR m_field.add_field("FIELD1", H1, AINSWORTH_LEGENDRE_BASE, 1);
252  CHKERR m_field.add_ents_to_field_by_type(one_prism_meshset, MBPRISM,
253  "FIELD1", VERY_NOISY);
254 
255  CHKERR m_field.set_field_order(one_prism_meshset, MBVERTEX, "FIELD1", 1);
256  CHKERR m_field.set_field_order(one_prism_meshset, MBEDGE, "FIELD1", 5,
257  VERY_NOISY);
258  CHKERR m_field.set_field_order(one_prism_meshset, MBTRI, "FIELD1", 5);
259  CHKERR m_field.set_field_order(one_prism_meshset, MBQUAD, "FIELD1", 5,
260  VERY_NOISY);
261  CHKERR m_field.set_field_order(one_prism_meshset, MBPRISM, "FIELD1", 7,
262  VERY_NOISY);
263  CHKERR m_field.build_fields(VERY_NOISY);
264 
265  // FE
266  CHKERR m_field.add_finite_element("PRISM");
267  CHKERR m_field.add_finite_element("EDGE");
268  CHKERR m_field.add_finite_element("TRI");
269  CHKERR m_field.add_finite_element("QUAD");
270 
271  // Define rows/cols and element data
272  CHKERR m_field.modify_finite_element_add_field_row("PRISM", "FIELD1");
273  CHKERR m_field.modify_finite_element_add_field_col("PRISM", "FIELD1");
274  CHKERR m_field.modify_finite_element_add_field_data("PRISM", "FIELD1");
275  // Define rows/cols and element data
276  CHKERR m_field.modify_finite_element_add_field_row("EDGE", "FIELD1");
277  CHKERR m_field.modify_finite_element_add_field_col("EDGE", "FIELD1");
278  CHKERR m_field.modify_finite_element_add_field_data("EDGE", "FIELD1");
279  // Define rows/cols and element data
280  CHKERR m_field.modify_finite_element_add_field_row("TRI", "FIELD1");
281  CHKERR m_field.modify_finite_element_add_field_col("TRI", "FIELD1");
282  CHKERR m_field.modify_finite_element_add_field_data("TRI", "FIELD1");
283  // Define rows/cols and element data
284  CHKERR m_field.modify_finite_element_add_field_row("QUAD", "FIELD1");
285  CHKERR m_field.modify_finite_element_add_field_col("QUAD", "FIELD1");
286  CHKERR m_field.modify_finite_element_add_field_data("QUAD", "FIELD1");
287 
288  CHKERR m_field.add_ents_to_finite_element_by_type(one_prism_meshset, MBEDGE,
289  "EDGE");
290  CHKERR m_field.add_ents_to_finite_element_by_type(one_prism_meshset, MBTRI,
291  "TRI");
292  CHKERR m_field.add_ents_to_finite_element_by_type(one_prism_meshset, MBQUAD,
293  "QUAD");
294  CHKERR m_field.add_ents_to_finite_element_by_type(one_prism_meshset,
295  MBPRISM, "PRISM");
296 
297  // build finite elemnts
298  CHKERR m_field.build_finite_elements();
299  // //build adjacencies
300  CHKERR m_field.build_adjacencies(bit_level0);
301 
302  // Problem
303  CHKERR m_field.add_problem("TEST_PROBLEM");
304 
305  // set finite elements for problem
306  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "PRISM");
307  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "EDGE");
308  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "TRI");
309  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "QUAD");
310  // set refinement level for problem
311  CHKERR m_field.modify_problem_ref_level_add_bit("TEST_PROBLEM", bit_level0);
312 
313  // build problem
314  ProblemsManager *prb_mng_ptr;
315  CHKERR m_field.getInterface(prb_mng_ptr);
316  CHKERR prb_mng_ptr->buildProblem("TEST_PROBLEM", true);
317  // partition
318  CHKERR prb_mng_ptr->partitionSimpleProblem("TEST_PROBLEM");
319  CHKERR prb_mng_ptr->partitionFiniteElements("TEST_PROBLEM");
320  // what are ghost nodes, see Petsc Manual
321  CHKERR prb_mng_ptr->partitionGhostDofs("TEST_PROBLEM");
322 
323  PrismFE fe_prism(m_field, tri_coords);
324  fe_prism.getOpPtrVector().push_back(new PrismOp(moab, map_coords));
325  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "PRISM", fe_prism);
326 
327  EdgeFE fe_edge(m_field, edge_block, one_prism);
328  fe_edge.getOpPtrVector().push_back(
330  moab, map_coords, one_prism));
331  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "EDGE", fe_edge);
332 
333  TriFE fe_tri(m_field, tri_coords, one_prism);
334  fe_tri.getOpPtrVector().push_back(
336  moab, map_coords, one_prism));
337  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "TRI", fe_tri);
338 
339  QuadFE fe_quad(m_field, edge_block, one_prism);
340  fe_quad.getOpPtrVector().push_back(
342  moab, map_coords, one_prism));
343  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "QUAD", fe_quad);
344 
345  CHKERR moab.write_file("prism_mesh.vtk", "VTK", "", &meshset, 1);
346  CHKERR moab.write_file("one_prism_mesh.vtk", "VTK", "", &one_prism_meshset,
347  1);
348  }
349  CATCH_ERRORS;
350 
352  return 0;
353 }
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
Problem manager is used to build and partition problems.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
static constexpr std::array< double, 3 > d4
multi_index_container< CoordsAndHandle, indexed_by< hashed_unique< composite_key< CoordsAndHandle, member< CoordsAndHandle, int, &CoordsAndHandle::x >, member< CoordsAndHandle, int, &CoordsAndHandle::y >, member< CoordsAndHandle, int, &CoordsAndHandle::z > > > > > MapCoords
virtual moab::Interface & get_moab()=0
static constexpr std::array< double, 3 > d3
MoFEMErrorCode setThickness(const Range &prisms, const double director3[], const double director4[])
MoFEMErrorCode getEntitiesByDimension(const int ms_id, const unsigned int cubit_bc_type, const int dimension, Range &entities, const bool recursive=true) const
get entities from CUBIT/meshset of a particular entity dimensionNodeset can contain nodes,...
MoFEMErrorCode seedPrismsEntities(Range &prisms, const BitRefLevel &bit, int verb=-1)
Seed prism entities by bit level.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string &name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
Interface for managing meshsets containing materials and boundary conditions.
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
Core (interface) class.
Definition: Core.hpp:50
char mesh_file_name[255]
static constexpr int number_of_prisms_layers
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:51
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
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
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
Managing BitRefLevels.
MoFEMErrorCode createPrisms(const Range &ents, Range &prisms, int verb=-1)
Make prisms from triangles.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:146
std::map< EntityHandle, EntityHandle > createdVertices
MoFEMErrorCode createPrismsFromPrisms(const Range &prisms, bool from_down, Range &out_prisms, int verb=-1)
Make prisms by extruding top or bottom prisms.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
static char help[]
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
#define CHKERR
Inline error check.
Definition: definitions.h:596
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:285
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elementsFunction which partition finite elements based on dofs partitioning....
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:433
continuous field
Definition: definitions.h:171
MoFEMErrorCode partitionSimpleProblem(const std::string name, int verb=VERBOSE)
partition problem dofs
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:61
DEPRECATED MoFEMErrorCode loop_finite_elements(const Problem *problem_ptr, const std::string &fe_name, FEMethod &method, int lower_rank, int upper_rank, MoFEMTypes bh, int verb=DEFAULT_VERBOSITY)

Variable Documentation

◆ d3

constexpr std::array<double, 3> d3 = {0, 0, 0}
static

◆ d4

constexpr std::array<double, 3> d4 = {0, 0, delta}
static

◆ delta

constexpr double delta
static
Initial value:
=
1. / static_cast<double>(number_of_prisms_layers)
static constexpr int number_of_prisms_layers
Examples
prism_elements_from_surface.cpp.

Definition at line 30 of file prism_elements_from_surface.cpp.

◆ help

char help[] = "...\n\n"
static

◆ number_of_prisms_layers

constexpr int number_of_prisms_layers = 18
static

◆ precision_exponent

constexpr int precision_exponent = 5
static