v0.14.0
Loading...
Searching...
No Matches
hcurl_check_approx_in_2d.cpp
Go to the documentation of this file.
1/**
2 * \file hcurl_check_approx_in_2d
3 * \example hcurl_check_approx_in_2d.cpp
4 *
5 * Approximate polynomial in 2D and check derivatives
6 *
7 */
8
9
10
11#include <MoFEM.hpp>
12
13using namespace MoFEM;
14
15static char help[] = "...\n\n";
16
19
20constexpr int BASE_DIM = 3;
21constexpr int SPACE_DIM = 2;
22
23/**
24 * @brief OPerator to integrate mass matrix for least square approximation
25 *
26 */
29
30/**
31 * @brief Operator to integrate the right hand side matrix for the problem
32 *
33 */
35 GAUSS>::OpSource<BASE_DIM, BASE_DIM>;
36
37constexpr double a0 = 0.0;
38constexpr double a1 = 2.0;
39constexpr double a2 = -15.0 * a0;
40constexpr double a3 = -20.0 / 6 * a1;
41constexpr double a4 = 15.0 * a0;
42constexpr double a5 = a1;
43constexpr double a6 = -a0;
45
46 static FTensor::Tensor1<double, BASE_DIM> fUn(const double x, const double y,
47 double z) {
49 // x
50 6 * a6 * std::pow(x, 5) * std::pow(y, 0) +
51 5 * a5 * std::pow(x, 4) * std::pow(y, 1) +
52 4 * a4 * std::pow(x, 3) * std::pow(y, 2) +
53 3 * a3 * std::pow(x, 2) * std::pow(y, 3) +
54 2 * a2 * std::pow(x, 1) * std::pow(y, 4) +
55 1 * a1 * std::pow(x, 0) * std::pow(y, 5),
56 // y
57 1 * a5 * std::pow(x, 5) * std::pow(y, 0) +
58 2 * a4 * std::pow(x, 4) * std::pow(y, 1) +
59 3 * a3 * std::pow(x, 3) * std::pow(y, 2) +
60 4 * a2 * std::pow(x, 2) * std::pow(y, 3) +
61 5 * a1 * std::pow(x, 1) * std::pow(y, 4) +
62 6 * a0 * std::pow(x, 0) * std::pow(y, 5),
63
64 // z
65 0.);
66 }
67
69 const double y) {
71 // x,x
72 30 * a6 * pow(x, 4) * pow(y, 0) + 20 * a5 * pow(x, 3) * pow(y, 1) +
73 12 * a4 * pow(x, 2) * pow(y, 2) + 6 * a3 * pow(x, 1) * pow(y, 3) +
74 2 * a2 * pow(x, 0) * pow(y, 4),
75 // x,y
76 5 * a5 * pow(x, 4) * pow(y, 0) + 8 * a4 * pow(x, 3) * pow(y, 1) +
77 9 * a3 * pow(x, 2) * pow(y, 2) + 8 * a2 * pow(x, 1) * pow(y, 3) +
78 5 * a1 * pow(x, 0) * pow(y, 4),
79 // y,x
80 5 * a5 * pow(x, 4) * pow(y, 0) + 8 * a4 * pow(x, 3) * pow(y, 1) +
81 9 * a3 * pow(x, 2) * pow(y, 2) + 8 * a2 * pow(x, 1) * pow(y, 3) +
82 5 * a1 * pow(x, 0) * pow(y, 4),
83 // y,y
84 2 * a4 * pow(x, 4) * pow(y, 0) + 6 * a3 * pow(x, 3) * pow(y, 1) +
85 12 * a2 * pow(x, 2) * pow(y, 2) + 20 * a1 * pow(x, 1) * pow(y, 3) +
86 30 * a0 * pow(x, 0) * pow(y, 4),
87 // z
88 0., 0.);
89 }
90
92 diff2Fun(const double x, const double y) {
94 // x,xx 0/000
95
96 30 * 4 * a6 * pow(x, 3) * pow(y, 0) +
97 20 * 3 * a5 * pow(x, 2) * pow(y, 1) +
98 12 * 2 * a4 * pow(x, 1) * pow(y, 2) +
99 6 * 1 * a3 * pow(x, 0) * pow(y, 3),
100
101 // x,xy 1/001
102
103 20 * 1 * a5 * pow(x, 3) * pow(y, 0) +
104 12 * 2 * a4 * pow(x, 2) * pow(y, 2) +
105 6 * 3 * a3 * pow(x, 1) * pow(y, 2) +
106 2 * 4 * a2 * pow(x, 0) * pow(y, 3),
107
108 // x,yx 2/010
109
110 5 * 4 * a5 * pow(x, 3) * pow(y, 0) +
111 8 * 3 * a4 * pow(x, 2) * pow(y, 1) +
112 9 * 2 * a3 * pow(x, 1) * pow(y, 2) +
113 8 * 1 * a2 * pow(x, 0) * pow(y, 3),
114
115 // x,yy 3/011
116
117 8 * 1 * a4 * pow(x, 3) * pow(y, 0) +
118 9 * 2 * a3 * pow(x, 2) * pow(y, 1) +
119 8 * 3 * a2 * pow(x, 1) * pow(y, 2) +
120 5 * 4 * a1 * pow(x, 0) * pow(y, 3),
121
122 // y,xx 4/100
123
124 5 * 4 * a5 * pow(x, 3) * pow(y, 0) +
125 8 * 3 * a4 * pow(x, 2) * pow(y, 1) +
126 9 * 2 * a3 * pow(x, 1) * pow(y, 2) +
127 8 * 1 * a2 * pow(x, 0) * pow(y, 3),
128
129 // y,xy 5/101
130
131 8 * 1 * a4 * pow(x, 3) * pow(y, 0) +
132 9 * 2 * a3 * pow(x, 2) * pow(y, 1) +
133 8 * 3 * a2 * pow(x, 1) * pow(y, 2) +
134 5 * 4 * a1 * pow(x, 0) * pow(y, 3),
135
136 // y,yx 6/110
137
138 2 * 4 * a4 * pow(x, 3) * pow(y, 0) +
139 6 * 3 * a3 * pow(x, 2) * pow(y, 1) +
140 12 * 2 * a2 * pow(x, 1) * pow(y, 2) +
141 20 * 1 * a1 * pow(x, 0) * pow(y, 3),
142
143 // y,yy 7/111
144
145 6 * 1 * a3 * pow(x, 3) * pow(y, 0) +
146 12 * 2 * a2 * pow(x, 2) * pow(y, 1) +
147 20 * 3 * a1 * pow(x, 1) * pow(y, 2) +
148 30 * 4 * a0 * pow(x, 0) * pow(y, 3),
149
150 // z,xx 8/200
151 0.,
152 // z,xy 9/201
153 0.,
154 // z,yx 10/210
155 0.,
156 // z,yy 11/211
157 0.);
158 }
159};
160
162 boost::shared_ptr<MatrixDouble> ptrVals;
163 boost::shared_ptr<VectorDouble> ptrDiv;
164 boost::shared_ptr<MatrixDouble> ptrGrad;
165 boost::shared_ptr<MatrixDouble> ptrHess;
166
167 OpCheckValsDiffVals(boost::shared_ptr<MatrixDouble> ptr_vals,
168 boost::shared_ptr<VectorDouble> ptr_div,
169 boost::shared_ptr<MatrixDouble> ptr_grad,
170 boost::shared_ptr<MatrixDouble> ptr_hess)
171 : FaceEleOp("FIELD1", OPROW), ptrVals(ptr_vals), ptrDiv(ptr_div),
172 ptrGrad(ptr_grad), ptrHess(ptr_hess) {}
173
177
181 const double eps = 1e-6;
182 if (type == MBEDGE && side == 0) {
183 const int nb_gauss_pts = data.getN().size1();
184
185 auto t_vals_from_op = getFTensor1FromMat<3>(*ptrVals);
186 auto t_div_from_op = getFTensor0FromVec(*ptrDiv);
187 auto t_grad_from_op = getFTensor2FromMat<3, 2>(*ptrGrad);
188 auto t_hess_from_op = getFTensor3FromMat<3, 2, 2>(*ptrHess);
189
190 for (int gg = 0; gg != nb_gauss_pts; gg++) {
191 const double x = getCoordsAtGaussPts()(gg, 0);
192 const double y = getCoordsAtGaussPts()(gg, 1);
193
194 // Check approximation
196 delta_val(i) = t_vals_from_op(i) - ApproxFunctions::fUn(x, y, 0)(i);
197 double err_val = sqrt(delta_val(i) * delta_val(i));
198 if (err_val > eps)
199 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
200 "Wrong value %4.3e", err_val);
201
202 FTensor::Tensor2<double, 3, 2> delta_diff_val;
203 delta_diff_val(i, j) =
204 t_grad_from_op(i, j) - ApproxFunctions::diffFun(x, y)(i, j);
205 double err_diff_val = sqrt(delta_diff_val(i, j) * delta_diff_val(i, j));
206 if (err_diff_val > eps)
207 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
208 "Wrong derivative of value %4.3e", err_diff_val);
209
210 double div = t_grad_from_op(0, 0) + t_grad_from_op(1, 1);
211 double err_div = div - t_div_from_op;
212 if (err_div > eps)
213 SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
214 "Wrong divergence from operator %4.3e (%4.3e != %4.3e)",
215 err_div, div, t_div_from_op);
216
217 FTensor::Tensor3<double, 3, 2, 2> delta_diff2_val;
218 delta_diff2_val(i, j, k) =
219 t_hess_from_op(i, j, k) - ApproxFunctions::diff2Fun(x, y)(i, j, k);
220 double hess_diff_error =
221 sqrt(delta_diff2_val(i, j, k) * delta_diff2_val(i, j, k));
222 if (hess_diff_error > eps)
223 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
224 "Wrong hessian from operator %4.3e", hess_diff_error);
225
226 ++t_vals_from_op;
227 ++t_div_from_op;
228 ++t_grad_from_op;
229 ++t_hess_from_op;
230 }
231 }
233 }
234};
235
236int main(int argc, char *argv[]) {
237
238 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
239
240 try {
241
242 DMType dm_name = "DMMOFEM";
243 CHKERR DMRegister_MoFEM(dm_name);
244
245 moab::Core mb_instance;
246 moab::Interface &moab = mb_instance;
247
248 // Create MoFEM instance
249 MoFEM::Core core(moab);
250 MoFEM::Interface &m_field = core;
251
252 Simple *simple_interface = m_field.getInterface<Simple>();
253 PipelineManager *pipeline_mng = m_field.getInterface<PipelineManager>();
254 CHKERR simple_interface->getOptions();
255 CHKERR simple_interface->loadFile("", "rectangle_tri.h5m");
256
257 // Declare elements
258 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOP };
259 const char *list_bases[] = {"ainsworth", "demkowicz"};
260 PetscBool flg;
261 PetscInt choice_base_value = AINSWORTH;
262 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
263 LASBASETOP, &choice_base_value, &flg);
264
265 if (flg != PETSC_TRUE)
266 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "base not set");
268 if (choice_base_value == AINSWORTH)
270 else if (choice_base_value == DEMKOWICZ)
272
273 CHKERR simple_interface->addDomainField("FIELD1", HCURL, base, 1);
274 constexpr int order = 5;
275 CHKERR simple_interface->setFieldOrder("FIELD1", order);
276 CHKERR simple_interface->setUp();
277 auto dm = simple_interface->getDM();
278
279 MatrixDouble vals, diff_vals;
280
281 auto assemble_matrices_and_vectors = [&]() {
283 auto jac_ptr = boost::make_shared<MatrixDouble>();
284 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
285 auto det_ptr = boost::make_shared<VectorDouble>();
286
287 pipeline_mng->getOpDomainRhsPipeline().push_back(
288 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
289 pipeline_mng->getOpDomainRhsPipeline().push_back(
290 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
291 pipeline_mng->getOpDomainRhsPipeline().push_back(
292 new OpMakeHdivFromHcurl());
293 pipeline_mng->getOpDomainRhsPipeline().push_back(
295 pipeline_mng->getOpDomainRhsPipeline().push_back(
296 new OpDomainSource("FIELD1", ApproxFunctions::fUn));
297
298 pipeline_mng->getOpDomainLhsPipeline().push_back(
299 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
300 pipeline_mng->getOpDomainLhsPipeline().push_back(
301 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
302 pipeline_mng->getOpDomainLhsPipeline().push_back(
303 new OpMakeHdivFromHcurl());
304 pipeline_mng->getOpDomainLhsPipeline().push_back(
306 pipeline_mng->getOpDomainLhsPipeline().push_back(
307
308 new OpDomainMass("FIELD1", "FIELD1",
309 [](double, double, double) { return 1; })
310
311 );
312
313 auto integration_rule = [](int, int, int p_data) { return 2 * p_data; };
314 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
315 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
316
318 };
319
320 auto solve_problem = [&] {
322 auto solver = pipeline_mng->createKSP();
323 CHKERR KSPSetFromOptions(solver);
324 CHKERR KSPSetUp(solver);
325
326 auto D = createDMVector(dm);
327 auto F = vectorDuplicate(D);
328
329 CHKERR KSPSolve(solver, F, D);
330 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
331 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
332 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
334 };
335
336 auto check_solution = [&] {
338
339 pipeline_mng->getOpDomainLhsPipeline().clear();
340 pipeline_mng->getOpDomainRhsPipeline().clear();
341
342 auto ptr_values = boost::make_shared<MatrixDouble>();
343 auto ptr_divergence = boost::make_shared<VectorDouble>();
344 auto ptr_grad = boost::make_shared<MatrixDouble>();
345 auto ptr_hessian = boost::make_shared<MatrixDouble>();
346
347 auto jac_ptr = boost::make_shared<MatrixDouble>();
348 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
349 auto det_ptr = boost::make_shared<VectorDouble>();
350
351 // Change H-curl to H-div in 2D, and apply Piola transform
352 pipeline_mng->getOpDomainRhsPipeline().push_back(
353 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
354 pipeline_mng->getOpDomainRhsPipeline().push_back(
355 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
356 pipeline_mng->getOpDomainRhsPipeline().push_back(
357 new OpMakeHdivFromHcurl());
358 pipeline_mng->getOpDomainRhsPipeline().push_back(
360 pipeline_mng->getOpDomainRhsPipeline().push_back(
361 new OpSetInvJacHcurlFace(inv_jac_ptr));
362
363 // Evaluate base function second derivative
364 auto base_mass = boost::make_shared<MatrixDouble>();
365 auto data_l2 = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
366
367 pipeline_mng->getOpDomainRhsPipeline().push_back(
368 new OpBaseDerivativesMass<BASE_DIM>(base_mass, data_l2, base, L2));
369 pipeline_mng->getOpDomainRhsPipeline().push_back(
371 inv_jac_ptr));
372 pipeline_mng->getOpDomainRhsPipeline().push_back(
373 new OpBaseDerivativesNext<BASE_DIM>(BaseDerivatives::SecondDerivative,
374 base_mass, data_l2, base, HCURL));
375
376 // Calculate field values at integration points
377 pipeline_mng->getOpDomainRhsPipeline().push_back(
378 new OpCalculateHVecVectorField<BASE_DIM>("FIELD1", ptr_values));
379 // Gradient
380 pipeline_mng->getOpDomainRhsPipeline().push_back(
382 ptr_grad));
383 // Hessian
384 pipeline_mng->getOpDomainRhsPipeline().push_back(
386 "FIELD1", ptr_divergence));
387 pipeline_mng->getOpDomainRhsPipeline().push_back(
389 ptr_hessian));
390
391 pipeline_mng->getOpDomainRhsPipeline().push_back(new OpCheckValsDiffVals(
392 ptr_values, ptr_divergence, ptr_grad, ptr_hessian));
393
394 CHKERR pipeline_mng->loopFiniteElements();
395
397 };
398
399 CHKERR assemble_matrices_and_vectors();
400 CHKERR solve_problem();
401 CHKERR check_solution();
402 }
404
406}
int main()
Definition: adol-c_atom.cpp:46
static const double eps
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ 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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
@ F
auto integration_rule
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:509
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
constexpr double a5
constexpr double a2
FormsIntegrators< FaceEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, BASE_DIM > OpDomainSource
Operator to integrate the right hand side matrix for the problem.
static char help[]
constexpr double a4
constexpr int SPACE_DIM
constexpr double a0
constexpr int BASE_DIM
constexpr double a3
constexpr double a1
constexpr double a6
FormsIntegrators< FaceEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, BASE_DIM > OpDomainMass
OPerator to integrate mass matrix for least square approximation.
double D
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
static FTensor::Tensor2< double, BASE_DIM, SPACE_DIM > diffFun(const double x, const double y)
static FTensor::Tensor3< double, BASE_DIM, SPACE_DIM, SPACE_DIM > diff2Fun(const double x, const double y)
static FTensor::Tensor1< double, BASE_DIM > fUn(const double x, const double y, double z)
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)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
@ OPROW
operator doWork function is executed on FE rows
Get vector field for H-div approximation.
Calculate gradient of vector field.
Calculate gradient of vector field.
Calculate divergence of vector field.
Make Hdiv space from Hcurl space in 2d.
PipelineManager interface.
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 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 getDM(DM *dm)
Get DM.
Definition: Simple.cpp:667
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
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.
FTensor::Index< 'i', 3 > i
boost::shared_ptr< MatrixDouble > ptrHess
boost::shared_ptr< MatrixDouble > ptrGrad
boost::shared_ptr< VectorDouble > ptrDiv
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCheckValsDiffVals(boost::shared_ptr< MatrixDouble > ptr_vals, boost::shared_ptr< VectorDouble > ptr_div, boost::shared_ptr< MatrixDouble > ptr_grad, boost::shared_ptr< MatrixDouble > ptr_hess)
FTensor::Index< 'j', 2 > j
boost::shared_ptr< MatrixDouble > ptrVals
FTensor::Index< 'k', 2 > k