v0.13.0
hcurl_divergence_operator_2d.cpp
Go to the documentation of this file.
1 /**
2  * \file hcurl_divergence_operator_2d.cpp
3  * \example hcurl_divergence_operator_2d.cpp
4  *
5  * Testing Hcurl base, transfromed to Hdiv base in 2d using Green theorem.
6  *
7  * Note 0: This is low-level implementation.
8  * Note 1: Generic implementation for Quad/Tri mesh of arbitrary order.
9  *
10  */
11 
12 /* This file is part of MoFEM.
13  * MoFEM is free software: you can redistribute it and/or modify it under
14  * the terms of the GNU Lesser General Public License as published by the
15  * Free Software Foundation, either version 3 of the License, or (at your
16  * option) any later version.
17  *
18  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
19  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
21  * License for more details.
22  *
23  * You should have received a copy of the GNU Lesser General Public
24  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
25 
26 #include <MoFEM.hpp>
27 
28 using namespace MoFEM;
29 
30 static char help[] = "...\n\n";
31 
34 
37 
38 constexpr int SPACE_DIM = 2;
40 
41 struct OpDivergence : public FaceEleOp {
42 
43  double &dIv;
44  OpDivergence(double &div) : FaceEleOp("FIELD1", OPROW), dIv(div) {}
45 
46  MoFEMErrorCode doWork(int side, EntityType type,
49  const int nb_dofs = data.getIndices().size();
50  if (nb_dofs == 0)
52  const int nb_gauss_pts = data.getN().size1();
53  auto t_w = getFTensor0IntegrationWeight();
54  auto t_diff_base_fun = data.getFTensor2DiffN<3, 2>();
55  const auto area = getMeasure();
56  for (int gg = 0; gg != nb_gauss_pts; gg++) {
57  const double val = area * t_w;
58  for (int bb = 0; bb != nb_dofs; bb++) {
59  dIv += val * t_diff_base_fun(i, i);
60  ++t_diff_base_fun;
61  }
62  ++t_w;
63  }
65  }
66 };
67 
68 struct OpFlux : public EdgeEleOp {
69 
70  double &fLux;
71  OpFlux(double &flux) : EdgeEleOp("FIELD1", OPROW), fLux(flux) {}
72 
73  MoFEMErrorCode doWork(int side, EntityType type,
76  const int nb_dofs = data.getIndices().size();
77  if (nb_dofs == 0)
79  const int nb_gauss_pts = data.getN().size1();
80  auto t_normal = getFTensor1Normal();
81  auto t_base_fun = data.getFTensor1N<3>();
82  auto t_w = getFTensor0IntegrationWeight();
83  for (int gg = 0; gg != nb_gauss_pts; gg++) {
84  for (int bb = 0; bb != nb_dofs; bb++) {
85  fLux += t_w * t_normal(i) * t_base_fun(i);
86  ++t_base_fun;
87  }
88  ++t_w;
89  }
90 
92  }
93 };
94 
95 int main(int argc, char *argv[]) {
96 
97  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
98 
99  try {
100 
101  moab::Core mb_instance;
102  moab::Interface &moab = mb_instance;
103 
104  PetscBool flg_file = PETSC_TRUE;
105  char mesh_file_name[255];
106  CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
107  255, &flg_file);
108  if (flg_file != PETSC_TRUE)
109  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
110  "*** ERROR -my_file (MESH FILE NEEDED)");
111 
112  // Read mesh to MOAB
113  CHKERR moab.load_file(mesh_file_name, 0, "");
114 
115  // Create MoFEM instance
116  MoFEM::Core core(moab);
117  MoFEM::Interface &m_field = core;
118 
119  // set entities bit level
120  BitRefLevel bit_level0 = BitRefLevel().set(0);
121  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
122  0, 2, bit_level0);
123 
124  // Declare elements
125  enum bases { AINSWORTH, DEMKOWICZ, LASBASETOP };
126  const char *list_bases[] = {"ainsworth", "demkowicz"};
127  PetscBool flg;
128  PetscInt choice_base_value = AINSWORTH;
129  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
130  LASBASETOP, &choice_base_value, &flg);
131  if (flg != PETSC_TRUE)
132  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "base not set");
134  if (choice_base_value == AINSWORTH)
136  else if (choice_base_value == DEMKOWICZ)
137  base = DEMKOWICZ_JACOBI_BASE;
138  int order = 5;
139  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
140 
141  CHKERR m_field.add_field("FIELD1", HCURL, base, 1);
142  CHKERR m_field.add_finite_element("FACE_FE");
143  // Define rows/cols and element data
144  CHKERR m_field.modify_finite_element_add_field_row("FACE_FE", "FIELD1");
145  CHKERR m_field.modify_finite_element_add_field_col("FACE_FE", "FIELD1");
146  CHKERR m_field.modify_finite_element_add_field_data("FACE_FE", "FIELD1");
147 
148  CHKERR m_field.add_finite_element("EDGE_FE");
149  // Define rows/cols and element data
150  CHKERR m_field.modify_finite_element_add_field_row("EDGE_FE", "FIELD1");
151  CHKERR m_field.modify_finite_element_add_field_col("EDGE_FE", "FIELD1");
152  CHKERR m_field.modify_finite_element_add_field_data("EDGE_FE", "FIELD1");
153 
154  // Problem
155  CHKERR m_field.add_problem("TEST_PROBLEM");
156  // set finite elements for problem
157  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "FACE_FE");
158  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "EDGE_FE");
159  // set refinement level for problem
160  CHKERR m_field.modify_problem_ref_level_add_bit("TEST_PROBLEM", bit_level0);
161 
162  // Add entities
163 
164  CHKERR m_field.add_ents_to_field_by_dim(0, SPACE_DIM, "FIELD1");
165  // Set order
166  for (auto t : {MBEDGE, MBTRI, MBQUAD})
167  CHKERR m_field.set_field_order(0, t, "FIELD1", order);
168 
169  // Add entities to elements
170  CHKERR m_field.add_ents_to_finite_element_by_dim(0, SPACE_DIM, "FACE_FE");
171  auto set_edge_elements_entities_on_mesh_skin = [&]() {
173  Range faces;
174  CHKERR moab.get_entities_by_dimension(0, 2, faces, false);
175  Skinner skin(&m_field.get_moab());
176  Range faces_skin;
177  CHKERR skin.find_skin(0, faces, false, faces_skin);
179  faces_skin, SPACE_DIM - 1, "EDGE_FE");
181  };
182  CHKERR set_edge_elements_entities_on_mesh_skin();
183 
184  // Build database
185  CHKERR m_field.build_fields();
186  // build finite elemnts
187  CHKERR m_field.build_finite_elements();
188  // build adjacencies
189  CHKERR m_field.build_adjacencies(bit_level0);
190 
191  // build problem
192  ProblemsManager *prb_mng_ptr;
193  CHKERR m_field.getInterface(prb_mng_ptr);
194  CHKERR prb_mng_ptr->buildProblem("TEST_PROBLEM", true);
195  // Partition
196  CHKERR prb_mng_ptr->partitionSimpleProblem("TEST_PROBLEM");
197  CHKERR prb_mng_ptr->partitionFiniteElements("TEST_PROBLEM");
198  CHKERR prb_mng_ptr->partitionGhostDofs("TEST_PROBLEM");
199 
200  // integration rule
201  auto rule = [&](int, int, int p) { return 2 * p; };
202 
203  auto calculate_divergence = [&]() {
204  double div = 0;
205  FaceEle fe_face(m_field);
206  fe_face.getRuleHook = rule;
207  auto jac_ptr = boost::make_shared<MatrixDouble>();
208  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
209  auto det_ptr = boost::make_shared<VectorDouble>();
210  fe_face.getOpPtrVector().push_back(new OpCalculateHOJacForFace(jac_ptr));
211  fe_face.getOpPtrVector().push_back(
212  new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
213  fe_face.getOpPtrVector().push_back(new OpMakeHdivFromHcurl());
214  fe_face.getOpPtrVector().push_back(
216  fe_face.getOpPtrVector().push_back(new OpSetInvJacHcurlFace(inv_jac_ptr));
217  fe_face.getOpPtrVector().push_back(
218  new OpSetHOWeightsOnFace());
219  fe_face.getOpPtrVector().push_back(new OpDivergence(div));
220  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "FACE_FE", fe_face);
221  return div;
222  };
223 
224  auto calculate_flux = [&]() {
225  double flux = 0;
226  EdgeEle fe_edge(m_field);
227  fe_edge.getRuleHook = rule;
228  fe_edge.getOpPtrVector().push_back(
230  fe_edge.getOpPtrVector().push_back(new OpFlux(flux));
231  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "EDGE_FE", fe_edge);
232  return flux;
233  };
234 
235  const double div = calculate_divergence();
236  const double flux = calculate_flux();
237 
238  MOFEM_LOG_CHANNEL("WOLD"); // reset channel
239  MOFEM_LOG_C("WORLD", Sev::inform,
240  "Div = %4.3e Flux = %3.4e Error = %4.3e\n", div, flux,
241  div - flux);
242 
243  constexpr double tol = 1e-8;
244  if (std::abs(div - flux) > tol)
245  SETERRQ2(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
246  "Test failed (div != flux) %3.4e != %3.4e", div, flux);
247  }
248  CATCH_ERRORS;
249 
251 }
static Index< 'p', 3 > p
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:314
MoFEM::FaceElementForcesAndSourcesCore FaceEle
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
FieldApproximationBase
approximation base
Definition: definitions.h:71
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:79
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ HCURL
field with continuous tangents
Definition: definitions.h:99
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_IMPOSIBLE_CASE
Definition: definitions.h:48
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:53
@ MOFEM_INVALID_DATA
Definition: definitions.h:49
#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
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
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
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 build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode add_ents_to_field_by_dim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
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.
EdgeElementForcesAndSourcesCoreSwitch< 0 > EdgeElementForcesAndSourcesCore
Edge finite element default.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:287
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
MoFEMErrorCode partitionSimpleProblem(const std::string name, int verb=VERBOSE)
partition problem dofs
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elements
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
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
int main(int argc, char *argv[])
static char help[]
FTensor::Index< 'i', SPACE_DIM > i
constexpr int SPACE_DIM
double tol
char mesh_file_name[255]
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
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
CoreTmp< 0 > Core
Definition: Core.hpp:1096
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
double flux
impulse magnitude
constexpr double t
plate stiffness
Definition: plate.cpp:76
Managing BitRefLevels.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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 moab::Interface & get_moab()=0
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
RuleHookFun getRuleHook
Hook to get rule.
Make Hdiv space from Hcurl space in 2d.
Modify integration weights on face to take in account higher-order geometry.
Problem manager is used to build and partition problems.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)