v0.14.0
Loading...
Searching...
No Matches
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
17using namespace MoFEM;
18
19static 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
45int 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
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 }
185
187}
188
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}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
int main()
Definition: adol-c_atom.cpp:46
static const double eps
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_IMPOSSIBLE_CASE
Definition: definitions.h:35
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr int order
auto integration_rule
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
virtual bool check_field(const std::string &name) const =0
check if field is in database
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.
static char help[]
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Core (interface) class.
Definition: Core.hpp:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
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 & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Calculate normals at Gauss points of triangle element.
transform Hcurl base fluxes from reference element to physical triangle
Apply covariant (Piola) transfer to Hcurl space for HO geometry.
transform local reference derivatives of shape function to global derivatives if higher order geometr...
Set inverse jacobian to base functions.
PipelineManager interface.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
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:320
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
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
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:282
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
OpFacesRot(FTensor::Tensor1< double, 3 > &t_curl)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
FTensor::Tensor1< double, 3 > & tCurl
OpVolCurl(FTensor::Tensor1< double, 3 > &t_curl)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
FTensor::Tensor1< double, 3 > & tCurl