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

Go to the source code of this file.

Classes

struct  MyOp
 
struct  CallingOp
 
struct  MyOp2
 

Typedefs

typedef tee_device< std::ostream, std::ofstream > TeeDevice
 
typedef stream< TeeDeviceTeeStream
 

Functions

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

Variables

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

Typedef Documentation

◆ TeeDevice

typedef tee_device<std::ostream, std::ofstream> TeeDevice

◆ TeeStream

typedef stream<TeeDevice> TeeStream

Function Documentation

◆ main()

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

Definition at line 180 of file forces_and_sources_testing_contact_prism_element.cpp.

180  {
181 
182  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
183 
184  try {
185 
186  moab::Core mb_instance;
187  moab::Interface &moab = mb_instance;
188  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
189  if (pcomm == NULL)
190  pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
191 
192  PetscBool flg = PETSC_TRUE;
193  char mesh_file_name[255];
194 #if PETSC_VERSION_GE(3, 6, 4)
195  CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
196  255, &flg);
197 #else
198  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
199  mesh_file_name, 255, &flg);
200 #endif
201  if (flg != PETSC_TRUE)
202  SETERRQ(PETSC_COMM_SELF, 1, "error -my_file (mesh file not given");
203 
204  const char *option;
205  option = "";
206  CHKERR moab.load_file(mesh_file_name, 0, option);
207 
208  // Create MoFEM (Joseph) database
209  MoFEM::Core core(moab);
210  MoFEM::Interface &m_field = core;
211  PrismInterface *interface;
212  CHKERR m_field.getInterface(interface);
213 
214  // set entities bit level
215  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
216  0, 3, BitRefLevel().set(0));
217  std::vector<BitRefLevel> bit_levels;
218  bit_levels.push_back(BitRefLevel().set(0));
219 
220  int ll = 1;
221  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, SIDESET, cit)) {
222 
223  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Insert Interface %d\n",
224  cit->getMeshsetId());
225  EntityHandle cubit_meshset = cit->getMeshset();
226  {
227  // get tet enties form back bit_level
228  EntityHandle ref_level_meshset = 0;
229  CHKERR moab.create_meshset(MESHSET_SET, ref_level_meshset);
231  ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
232  BitRefLevel().set(), MBTET,
233  ref_level_meshset);
235  ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
236  BitRefLevel().set(), MBPRISM,
237  ref_level_meshset);
238  Range ref_level_tets;
239  CHKERR moab.get_entities_by_handle(ref_level_meshset, ref_level_tets,
240  true);
241  // get faces and test to split
242  CHKERR interface->getSides(cubit_meshset, bit_levels.back(), true, 0);
243  // set new bit level
244  bit_levels.push_back(BitRefLevel().set(ll++));
245  // split faces and
246  CHKERR interface->splitSides(ref_level_meshset, bit_levels.back(),
247  cubit_meshset, true, true, 0);
248  // clean meshsets
249  CHKERR moab.delete_entities(&ref_level_meshset, 1);
250  }
251  // update cubit meshsets
252  for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, ciit)) {
253  EntityHandle cubit_meshset = ciit->meshset;
255  ->updateMeshsetByEntitiesChildren(cubit_meshset, bit_levels.back(),
256  cubit_meshset, MBVERTEX, true);
258  ->updateMeshsetByEntitiesChildren(cubit_meshset, bit_levels.back(),
259  cubit_meshset, MBEDGE, true);
261  ->updateMeshsetByEntitiesChildren(cubit_meshset, bit_levels.back(),
262  cubit_meshset, MBTRI, true);
264  ->updateMeshsetByEntitiesChildren(cubit_meshset, bit_levels.back(),
265  cubit_meshset, MBTET, true);
266  }
267  }
268 
269  // Fields
270  CHKERR m_field.add_field("FIELD1", H1, AINSWORTH_LEGENDRE_BASE, 3);
271  CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
272  3);
273  CHKERR m_field.add_field("FIELD2", NOFIELD, NOBASE, 3);
274 
275  auto set_no_field_vertex = [&]() {
277  // Creating and adding no field entities.
278  const double coords[] = {0, 0, 0};
279  EntityHandle no_field_vertex;
280  CHKERR m_field.get_moab().create_vertex(coords, no_field_vertex);
281  Range range_no_field_vertex;
282  range_no_field_vertex.insert(no_field_vertex);
283  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
284  range_no_field_vertex, BitRefLevel().set());
285  EntityHandle meshset = m_field.get_field_meshset("FIELD2");
286  CHKERR m_field.get_moab().add_entities(meshset, range_no_field_vertex);
288  };
289 
290  CHKERR set_no_field_vertex();
291 
292  // FE
293  CHKERR m_field.add_finite_element("TEST_FE1");
294  CHKERR m_field.add_finite_element("TEST_FE2");
295 
296  // Define rows/cols and element data
297  CHKERR m_field.modify_finite_element_add_field_row("TEST_FE1", "FIELD1");
298  CHKERR m_field.modify_finite_element_add_field_col("TEST_FE1", "FIELD1");
299  CHKERR m_field.modify_finite_element_add_field_data("TEST_FE1", "FIELD1");
300  CHKERR m_field.modify_finite_element_add_field_data("TEST_FE1",
301  "MESH_NODE_POSITIONS");
302 
303  CHKERR m_field.modify_finite_element_add_field_row("TEST_FE2", "FIELD1");
304  CHKERR m_field.modify_finite_element_add_field_col("TEST_FE2", "FIELD2");
305  CHKERR m_field.modify_finite_element_add_field_data("TEST_FE2", "FIELD1");
306  CHKERR m_field.modify_finite_element_add_field_data("TEST_FE2", "FIELD2");
307 
308  // Problem
309  CHKERR m_field.add_problem("TEST_PROBLEM");
310 
311  // set finite elements for problem
312  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM",
313  "TEST_FE1");
314  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM",
315  "TEST_FE2");
316  // set refinement level for problem
317  CHKERR m_field.modify_problem_ref_level_add_bit("TEST_PROBLEM",
318  bit_levels.back());
319 
320  // meshset consisting all entities in mesh
321  EntityHandle root_set = moab.get_root_set();
322  // add entities to field
323  CHKERR m_field.add_ents_to_field_by_type(root_set, MBTET, "FIELD1");
324  CHKERR m_field.add_ents_to_field_by_type(root_set, MBTET,
325  "MESH_NODE_POSITIONS");
326  // add entities to finite element
327  CHKERR m_field.add_ents_to_finite_element_by_type(root_set, MBPRISM,
328  "TEST_FE1", 10);
329  CHKERR m_field.add_ents_to_finite_element_by_type(root_set, MBPRISM,
330  "TEST_FE2", 10);
331 
332  // set app. order
333  // see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes
334  // (Mark Ainsworth & Joe Coyle)
335  int order = 3;
336  CHKERR m_field.set_field_order(root_set, MBTET, "FIELD1", order);
337  CHKERR m_field.set_field_order(root_set, MBTRI, "FIELD1", order);
338  CHKERR m_field.set_field_order(root_set, MBEDGE, "FIELD1", order);
339  CHKERR m_field.set_field_order(root_set, MBVERTEX, "FIELD1", 1);
340 
341  CHKERR m_field.set_field_order(root_set, MBTET, "MESH_NODE_POSITIONS", 2);
342  CHKERR m_field.set_field_order(root_set, MBTRI, "MESH_NODE_POSITIONS", 2);
343  CHKERR m_field.set_field_order(root_set, MBEDGE, "MESH_NODE_POSITIONS", 2);
344  CHKERR m_field.set_field_order(root_set, MBVERTEX, "MESH_NODE_POSITIONS",
345  1);
346 
347  /****/
348  // build database
349  // build field
350  CHKERR m_field.build_fields();
351  // set FIELD1 from positions of 10 node tets
352  Projection10NodeCoordsOnField ent_method_field1(m_field, "FIELD1");
353  CHKERR m_field.loop_dofs("FIELD1", ent_method_field1);
354  Projection10NodeCoordsOnField ent_method_mesh_positions(
355  m_field, "MESH_NODE_POSITIONS");
356  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_mesh_positions);
357 
358  // build finite elemnts
359  CHKERR m_field.build_finite_elements();
360  // build adjacencies
361  CHKERR m_field.build_adjacencies(bit_levels.back());
362  // build problem
363  ProblemsManager *prb_mng_ptr;
364  CHKERR m_field.getInterface(prb_mng_ptr);
365  CHKERR prb_mng_ptr->buildProblem("TEST_PROBLEM", false);
366 
367  /****/
368  // mesh partitioning
369  // partition
370  CHKERR prb_mng_ptr->partitionSimpleProblem("TEST_PROBLEM");
371  CHKERR prb_mng_ptr->partitionFiniteElements("TEST_PROBLEM");
372  // what are ghost nodes, see Petsc Manual
373  CHKERR prb_mng_ptr->partitionGhostDofs("TEST_PROBLEM");
374 
375 
376  std::ofstream ofs("forces_and_sources_testing_contact_prism_element.txt");
377  TeeDevice my_tee(std::cout, ofs);
378  TeeStream my_split(my_tee);
379 
381  using ContactDataOp =
383 
385  fe1.getOpPtrVector().push_back(
386  new MyOp(my_split, UMDataOp::OPROW, ContactDataOp::FACEMASTER));
387  fe1.getOpPtrVector().push_back(
388  new MyOp(my_split, UMDataOp::OPROW, ContactDataOp::FACESLAVE));
389  fe1.getOpPtrVector().push_back(new MyOp(my_split, UMDataOp::OPROWCOL,
390  ContactDataOp::FACEMASTERMASTER));
391  fe1.getOpPtrVector().push_back(
392  new MyOp(my_split, UMDataOp::OPROWCOL, ContactDataOp::FACEMASTERSLAVE));
393  fe1.getOpPtrVector().push_back(
394  new MyOp(my_split, UMDataOp::OPROWCOL, ContactDataOp::FACESLAVEMASTER));
395  fe1.getOpPtrVector().push_back(
396  new MyOp(my_split, UMDataOp::OPROWCOL, ContactDataOp::FACESLAVESLAVE));
397  fe1.getOpPtrVector().push_back(new CallingOp(my_split, UMDataOp::OPCOL));
398  fe1.getOpPtrVector().push_back(new CallingOp(my_split, UMDataOp::OPROW));
399  fe1.getOpPtrVector().push_back(new CallingOp(my_split, UMDataOp::OPROWCOL));
400  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "TEST_FE1", fe1);
401 
403  fe2.getOpPtrVector().push_back(
404  new MyOp2(my_split, UMDataOp::OPCOL, ContactDataOp::FACEMASTER));
405  fe2.getOpPtrVector().push_back(
406  new MyOp2(my_split, UMDataOp::OPCOL, ContactDataOp::FACESLAVE));
407  fe2.getOpPtrVector().push_back(new MyOp2(my_split, UMDataOp::OPROWCOL,
408  ContactDataOp::FACEMASTERMASTER));
409  fe2.getOpPtrVector().push_back(new MyOp2(my_split, UMDataOp::OPROWCOL,
410  ContactDataOp::FACEMASTERSLAVE));
411  fe2.getOpPtrVector().push_back(new MyOp2(my_split, UMDataOp::OPROWCOL,
412  ContactDataOp::FACESLAVEMASTER));
413  fe2.getOpPtrVector().push_back(
414  new MyOp2(my_split, UMDataOp::OPROWCOL, ContactDataOp::FACESLAVESLAVE));
415  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "TEST_FE2", fe2);
416  }
417  CATCH_ERRORS;
418 
420 
421  return 0;
422 }
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
Deprecated interface functions.
Problem manager is used to build and partition problems.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual moab::Interface & get_moab()=0
scalar or vector of scalars describe (no true field)
Definition: definitions.h:175
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
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
MoFEMErrorCode splitSides(const EntityHandle meshset, const BitRefLevel &bit, const int msId, const CubitBCType cubit_bc_type, const bool add_interface_entities, const bool recursive=false, int verb=QUIET)
Split nodes and other entities of tetrahedra on both sides of the interface and insert flat prisms in...
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.
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]
Projection of edge entities with one mid-node on hierarchical basis.
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.
virtual EntityHandle get_field_meshset(const std::string &name) const =0
get field meshset
constexpr int order
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.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:151
Create interface from given surface and insert flat prisms in-between.
tee_device< std::ostream, std::ofstream > TeeDevice
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
MoFEMErrorCode getSides(const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
Store tetrahedra from each side of the interface separately in two child meshsets of the parent meshs...
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
#define CHKERR
Inline error check.
Definition: definitions.h:601
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
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:290
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
ContactPrism finite elementUser is implementing own operator at Gauss points level,...
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
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 MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:438
continuous field
Definition: definitions.h:176
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

◆ help

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