v0.14.0
boundary_marker.cpp
Go to the documentation of this file.
1 /**
2  * \file boundary_marker.cpp
3  * \example boundary_marker.cpp
4  *
5  * Mark skin entities, i.e. boundary, next check DOFs if are properly marked
6  * when accessed form user data operator.
7  *
8  */
9 
10 #include <MoFEM.hpp>
11 
12 using namespace MoFEM;
13 
14 static char help[] = "...\n\n";
15 
18 
19 struct OpFace : public FaceEleOp {
20  const Range &skinEnts;
21  const std::vector<unsigned char> &mArker;
22 
23  OpFace(const Range &skin_ents, const std::vector<unsigned char> &marker)
24  : FaceEleOp("FIELD1", OPROW), skinEnts(skin_ents), mArker(marker) {}
25 
26  MoFEMErrorCode doWork(int side, EntityType type,
29 
30  const int nb_dofs = data.getIndices().size();
31  if (nb_dofs == 0)
33 
34  // This function takes entities on DOFs, and check if those entities are
35  // contained in the range. Note it is local element vector, which is used
36  // to validate of global local vector.
37  auto get_boundary_marker_directly_from_range = [&]() {
38  std::vector<unsigned char> ents_marker_used_for_testing;
39  ents_marker_used_for_testing.resize(data.getLocalIndices().size(),0);
40  switch (type) {
41  case MBVERTEX: {
42  for (size_t side = 0; side != getNumberOfNodesOnElement(); ++side) {
43  EntityHandle ent = getSideEntity(side, MBVERTEX);
44  ents_marker_used_for_testing[3 * side + 0] =
45  skinEnts.find(ent) != skinEnts.end() ? 1 : 0;
46  ents_marker_used_for_testing[3 * side + 2] =
47  ents_marker_used_for_testing[3 * side + 0];
48  }
49  break;
50  }
51  default: {
52  EntityHandle ent = getSideEntity(side, type);
53  bool is_on_boundary = skinEnts.find(ent) != skinEnts.end();
54  for (size_t dof = 0; dof != nb_dofs; ++dof)
55  if ((dof % 3) != 1)
56  ents_marker_used_for_testing[dof] = is_on_boundary ? 1 : 0;
57  }
58  };
59  return ents_marker_used_for_testing;
60  };
61 
62  auto test_marker = get_boundary_marker_directly_from_range();
63  for (size_t n = 0; n != nb_dofs; ++n) {
64  const size_t local_index = data.getLocalIndices()[n];
65  if (test_marker[n] != mArker[local_index])
66  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
67  "Wrong boundary marker");
68  }
69 
71  }
72 };
73 
74 int main(int argc, char *argv[]) {
75 
76  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
77 
78  try {
79 
80  moab::Core mb_instance;
81  moab::Interface &moab = mb_instance;
82 
83  // Create MoFEM instance
84  MoFEM::Core core(moab);
85  MoFEM::Interface &m_field = core;
86 
87  int order = 4;
88  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
89 
90  DMType dm_name = "DMMOFEM";
91  CHKERR DMRegister_MoFEM(dm_name);
92 
93  auto *simple_interface = m_field.getInterface<Simple>();
94  CHKERR simple_interface->getOptions();
95  CHKERR simple_interface->loadFile("");
96  CHKERR simple_interface->addDomainField("FIELD1", H1,
98  CHKERR simple_interface->setFieldOrder("FIELD1", order);
99  CHKERR simple_interface->setUp();
100 
101  auto get_ents_on_mesh_skin = [&]() {
102  Range faces;
103  CHKERR m_field.get_moab().get_entities_by_type(0, MBTRI, faces);
104  Skinner skin(&m_field.get_moab());
105  Range skin_edges;
106  CHKERR skin.find_skin(0, faces, false, skin_edges);
107  Range skin_verts;
108  CHKERR moab.get_connectivity(skin_edges, skin_verts, true);
109  skin_edges.merge(skin_verts);
110  return skin_edges;
111  };
112 
113  auto mark_boundary_dofs = [&](Range &skin_edges) {
114  auto problem_manager = m_field.getInterface<ProblemsManager>();
115  std::vector<unsigned char> marker;
116  problem_manager->markDofs(simple_interface->getProblemName(), ROW,
117  ProblemsManager::OR, skin_edges, marker);
118  // Unset coefficient 1, e.g. u_y
119  problem_manager->modifyMarkDofs(simple_interface->getProblemName(), ROW,
120  "FIELD1", 1, 1, ProblemsManager::AND, 0,
121  marker);
122  return marker;
123  };
124 
125  auto skin_ents = get_ents_on_mesh_skin();
126 
127  // Get global local vector of marked DOFs. Is global, since is set for all
128  // DOFs on processor. Is local since only DOFs on processor are in the
129  // vector. To access DOFs use local indices.
130  auto marker = mark_boundary_dofs(skin_ents);
131 
132  boost::shared_ptr<FaceEle> fe(new FaceEle(m_field));
133  fe->getOpPtrVector().push_back(new OpFace(skin_ents, marker));
134 
135  auto dm = simple_interface->getDM();
136  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
137  fe);
138  }
139  CATCH_ERRORS;
140 
142 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
EntityHandle
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
MoFEM::EntitiesFieldData::EntData::getLocalIndices
const VectorInt & getLocalIndices() const
get local indices of dofs on entity
Definition: EntitiesFieldData.hpp:1229
OpFace::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: boundary_marker.cpp:26
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::ProblemsManager::OR
@ OR
Definition: ProblemsManager.hpp:416
MoFEM.hpp
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ROW
@ ROW
Definition: definitions.h:136
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: FaceElementForcesAndSourcesCore.hpp:86
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
OpFace::OpFace
OpFace(const Range &skin_ents, const std::vector< unsigned char > &marker)
Definition: boundary_marker.cpp:23
convert.type
type
Definition: convert.py:64
help
static char help[]
Definition: boundary_marker.cpp:14
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1214
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
OpFace
Definition: boundary_marker.cpp:19
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MoFEM::ProblemsManager::AND
@ AND
Definition: ProblemsManager.hpp:416
convert.n
n
Definition: convert.py:82
FaceEle
MoFEM::FaceElementForcesAndSourcesCore FaceEle
Definition: boundary_marker.cpp:16
Range
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
main
int main(int argc, char *argv[])
Definition: boundary_marker.cpp:74
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
OpFace::skinEnts
const Range & skinEnts
Definition: boundary_marker.cpp:20
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
OpFace::mArker
const std::vector< unsigned char > & mArker
Definition: boundary_marker.cpp:21
marker
auto marker
set bit to marker
Definition: hanging_node_approx.cpp:82
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359