v0.14.0
hcurl_curl_operator.cpp
Go to the documentation of this file.
1 /**
2  * \file hcurl_curl_operator.cpp
3  * \brief Testich curl-curl operator by applying Stokes theorem
4  * \example hcurl_curl_operator.cpp
5  *
6  * Using PipelineManager interface calculate the curl of base functions, and
7  * integral of the vector tangent vector with normal on the boundary. Since the
8  * h-curl space is used, volume integral and boundary integral should give the
9  * same result, as a result, as we are applying Stokes theorem on h-curl space.
10  *
11  */
12 
13 
14 
15 #include <MoFEM.hpp>
16 
17 using namespace MoFEM;
18 
19 static char help[] = "...\n\n";
20 
22 
26  "HCURL", UserDataOperator::OPROW),
27  tCurl(t_curl) {}
28 
29  MoFEMErrorCode doWork(int side, EntityType type,
31 };
32 
34 
38  "HCURL", UserDataOperator::OPROW),
39  tCurl(t_curl) {}
40 
41  MoFEMErrorCode doWork(int side, EntityType type,
43 };
44 
45 int main(int argc, char *argv[]) {
46 
47  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
48 
49  try {
50 
51  enum bases { AINSWORTH, DEMKOWICZ, LASTOP };
52 
53  const char *list[] = {"ainsworth", "demkowicz"};
54 
55  PetscBool flg;
56  PetscInt choise_value = AINSWORTH;
57  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list, LASTOP,
58  &choise_value, &flg);
59  if (flg != PETSC_TRUE) {
60  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "base not set");
61  }
62 
63  PetscBool ho_geometry = PETSC_FALSE;
64  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-ho_geometry", &ho_geometry,
65  PETSC_NULL);
66 
67  DMType dm_name = "DMMOFEM";
68  CHKERR DMRegister_MoFEM(dm_name);
69 
70  // Create MoAB
71  moab::Core mb_instance; ///< database
72  moab::Interface &moab = mb_instance; ///< interface
73 
74  // Create MoFEM
75  MoFEM::Core core(moab);
76  MoFEM::Interface &m_field = core;
77 
78  Simple *simple_interface = m_field.getInterface<Simple>();
79  PipelineManager *pipeline_mng = m_field.getInterface<PipelineManager>();
80  CHKERR simple_interface->getOptions();
81  CHKERR simple_interface->loadFile("");
82 
83  // fields
84  switch (choise_value) {
85  case AINSWORTH:
86  CHKERR simple_interface->addDomainField("HCURL", HCURL,
88  CHKERR simple_interface->addBoundaryField("HCURL", HCURL,
90  break;
91  case DEMKOWICZ:
92  CHKERR simple_interface->addDomainField("HCURL", HCURL,
94  CHKERR simple_interface->addBoundaryField("HCURL", HCURL,
96  break;
97  }
98 
99  if (ho_geometry == PETSC_TRUE)
100  CHKERR simple_interface->addDataField("MESH_NODE_POSITIONS", H1,
102 
103  constexpr int order = 3;
104  CHKERR simple_interface->setFieldOrder("HCURL", order);
105  // Range ents;
106  // CHKERR moab.get_entities_by_dimension(0, 2, ents, true);
107  // CHKERR simple_interface->setFieldOrder("HCURL", 1, &ents);
108 
109  if (ho_geometry == PETSC_TRUE)
110  CHKERR simple_interface->setFieldOrder("MESH_NODE_POSITIONS", 2);
111  CHKERR simple_interface->setUp();
112 
113  auto integration_rule = [](int, int, int p_data) { return 2 * p_data; };
114  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
115  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
116 
117  FTensor::Tensor1<double, 3> t_curl_vol;
118  FTensor::Tensor1<double, 3> t_curl_skin;
119 
120  auto material_grad_mat = boost::make_shared<MatrixDouble>();
121  auto material_det_vec = boost::make_shared<VectorDouble>();
122  auto material_inv_grad_mat = boost::make_shared<MatrixDouble>();
123  auto jac_ptr = boost::make_shared<MatrixDouble>();
124  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
125  auto det_ptr = boost::make_shared<VectorDouble>();
126 
127  pipeline_mng->getOpDomainRhsPipeline().push_back(
128  new OpCalculateHOJac<3>(jac_ptr));
129  pipeline_mng->getOpDomainRhsPipeline().push_back(
130  new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
131  pipeline_mng->getOpDomainRhsPipeline().push_back(
132  new OpSetHOWeights(det_ptr));
133  pipeline_mng->getOpDomainRhsPipeline().push_back(
134  new OpSetHOCovariantPiolaTransform(HCURL, inv_jac_ptr));
135  pipeline_mng->getOpDomainRhsPipeline().push_back(
136  new OpSetHOInvJacVectorBase(HCURL, inv_jac_ptr));
137 
138  if (ho_geometry) {
139  pipeline_mng->getOpDomainRhsPipeline().push_back(
140  new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
141  material_grad_mat));
142  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpInvertMatrix<3>(
143  material_grad_mat, material_det_vec, material_inv_grad_mat));
144  pipeline_mng->getOpDomainRhsPipeline().push_back(
145  new OpSetHOWeights(material_det_vec));
146  pipeline_mng->getOpDomainRhsPipeline().push_back(
147  new OpSetHOCovariantPiolaTransform(HCURL, material_inv_grad_mat));
148  pipeline_mng->getOpDomainRhsPipeline().push_back(
149  new OpSetHOInvJacVectorBase(HCURL, material_inv_grad_mat));
150  }
151  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpVolCurl(t_curl_vol));
152 
153  if (m_field.check_field("MESH_NODE_POSITIONS"))
154  pipeline_mng->getOpBoundaryRhsPipeline().push_back(
155  new OpGetHONormalsOnFace("MESH_NODE_POSITIONS"));
156  pipeline_mng->getOpBoundaryRhsPipeline().push_back(
158  pipeline_mng->getOpBoundaryRhsPipeline().push_back(
159  new OpFacesRot(t_curl_skin));
160 
162  t_curl_vol(i) = 0;
163  t_curl_skin(i) = 0;
164 
165  // project geometry form 10 node tets on higher order approx. functions
166  if (ho_geometry == PETSC_TRUE) {
167  Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
168  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
169  }
170 
171  // Run pipelines on mesh
172  CHKERR pipeline_mng->loopFiniteElements();
173 
174  std::cout.precision(12);
175 
176  t_curl_vol(i) -= t_curl_skin(i);
177  double nrm2 = sqrt(t_curl_vol(i) * t_curl_vol(i));
178 
179  constexpr double eps = 1e-8;
180  if (fabs(nrm2) > eps)
181  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
182  "Curl operator not passed test\n");
183  }
184  CATCH_ERRORS;
185 
187 }
188 
189 MoFEMErrorCode OpVolCurl::doWork(int side, EntityType type,
192 
193  if (data.getFieldData().size() == 0)
195 
196  const unsigned int nb_gauss_pts = data.getDiffN().size1();
197  const unsigned int nb_dofs = data.getFieldData().size();
198 
199  MatrixDouble curl_mat;
203 
204  auto t_curl_base = data.getFTensor2DiffN<3, 3>();
205 
206  unsigned int gg = 0;
207  for (; gg < nb_gauss_pts; gg++) {
208  double w = getGaussPts()(3, gg) * getVolume();
210  for (unsigned int dd = 0; dd != nb_dofs; dd++) {
211  t_curl(i) = levi_civita(j, i, k) * t_curl_base(j, k);
212  tCurl(i) += w * t_curl(i);
213  ++t_curl_base;
214  }
215  }
216 
218 }
219 
223 
224  int nb_dofs = data.getFieldData().size();
225  if (nb_dofs == 0)
227  int nb_gauss_pts = data.getN().size1();
228 
229  auto t_curl_base = data.getFTensor1N<3>();
230 
233 
234  for (int gg = 0; gg < nb_gauss_pts; gg++) {
235  for (int dd = 0; dd < nb_dofs; dd++) {
236  double w = getGaussPts()(2, gg);
237  const double n0 = getNormalsAtGaussPts(gg)[0];
238  const double n1 = getNormalsAtGaussPts(gg)[1];
239  const double n2 = getNormalsAtGaussPts(gg)[2];
240  if (getFEType() == MBTRI) {
241  w *= 0.5;
242  }
243 
244  tCurl(0) += (n1 * t_curl_base(2) - n2 * t_curl_base(1)) * w;
245  tCurl(1) += (n2 * t_curl_base(0) - n0 * t_curl_base(2)) * w;
246  tCurl(2) += (n0 * t_curl_base(1) - n1 * t_curl_base(0)) * w;
247  ++t_curl_base;
248  }
249  }
250 
252 }
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::CoreInterface::loop_dofs
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.
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
FTensor::Tensor1< double, 3 >
OpVolCurl::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: hcurl_curl_operator.cpp:189
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::EntitiesFieldData::EntData::getFieldData
const VectorDouble & getFieldData() const
get dofs values
Definition: EntitiesFieldData.hpp:1254
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
OpFacesRot::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: hcurl_curl_operator.cpp:220
MoFEM.hpp
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::EntitiesFieldData::EntData::getDiffN
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
Definition: EntitiesFieldData.hpp:1329
OpVolCurl::tCurl
FTensor::Tensor1< double, 3 > & tCurl
Definition: hcurl_curl_operator.cpp:23
help
static char help[]
Definition: hcurl_curl_operator.cpp:19
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition: Levi_Civita.hpp:617
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MOFEM_IMPOSSIBLE_CASE
@ MOFEM_IMPOSSIBLE_CASE
Definition: definitions.h:35
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEM::EntitiesFieldData::EntData::getFTensor2DiffN
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
Definition: EntitiesFieldData.cpp:660
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::OpSetHOCovariantPiolaTransform
Apply covariant (Piola) transfer to Hcurl space for HO geometry.
Definition: HODataOperators.hpp:206
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
OpFacesRot::tCurl
FTensor::Tensor1< double, 3 > & tCurl
Definition: hcurl_curl_operator.cpp:35
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
convert.type
type
Definition: convert.py:64
main
int main(int argc, char *argv[])
Definition: hcurl_curl_operator.cpp:45
MoFEM::Simple::addDomainField
MoFEMErrorCode addDomainField(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_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:264
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MoFEM::OpCalculateHOJac< 3 >
Definition: HODataOperators.hpp:269
OpFacesRot
Definition: hcurl_curl_operator.cpp:33
MoFEM::OpGetHONormalsOnFace
Calculate normals at Gauss points of triangle element.
Definition: HODataOperators.hpp:281
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::Simple::addDataField
MoFEMErrorCode addDataField(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_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:392
OpVolCurl::OpVolCurl
OpVolCurl(FTensor::Tensor1< double, 3 > &t_curl)
Definition: hcurl_curl_operator.cpp:24
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:108
FTensor::Index< 'i', 3 >
OpFacesRot::OpFacesRot
OpFacesRot(FTensor::Tensor1< double, 3 > &t_curl)
Definition: hcurl_curl_operator.cpp:36
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:545
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1535
integration_rule
auto integration_rule
Definition: free_surface.cpp:185
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
FTensor::dd
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
MoFEM::OpHOSetCovariantPiolaTransformOnFace3D
transform Hcurl base fluxes from reference element to physical triangle
Definition: HODataOperators.hpp:336
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: EntitiesFieldData.hpp:1318
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
MoFEM::EntitiesFieldData::EntData::getFTensor1N
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Definition: EntitiesFieldData.cpp:640
MoFEM::OpSetHOWeights
Set inverse jacobian to base functions.
Definition: HODataOperators.hpp:159
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:203
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3249
MoFEM::Simple::addBoundaryField
MoFEMErrorCode addBoundaryField(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_ZERO, int verb=-1)
Add field on boundary.
Definition: Simple.cpp:354
sdf_wavy_2d.w
int w
Definition: sdf_wavy_2d.py:6
MoFEM::OpSetHOInvJacVectorBase
transform local reference derivatives of shape function to global derivatives if higher order geometr...
Definition: HODataOperators.hpp:91
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
OpVolCurl
Definition: hcurl_curl_operator.cpp:21
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:683
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
FractureMechanics::LASTOP
@ LASTOP
Definition: CrackFrontElement.hpp:23
convert.int
int
Definition: convert.py:64
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
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182