v0.13.1
Classes | Typedefs | Functions | Variables
hcurl_check_approx_in_2d.cpp File Reference
#include <MoFEM.hpp>

Go to the source code of this file.

Classes

struct  ApproxFunctions
 
struct  OpCheckValsDiffVals
 

Typedefs

using FaceEle = MoFEM::FaceElementForcesAndSourcesCore
 
using FaceEleOp = FaceEle::UserDataOperator
 
using OpDomainMass = FormsIntegrators< FaceEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, BASE_DIM >
 OPerator to integrate mass matrix for least square approximation. More...
 
using OpDomainSource = FormsIntegrators< FaceEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, BASE_DIM >
 Operator to integrate the right hand side matrix for the problem. More...
 

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 
constexpr int BASE_DIM = 3
 
constexpr int SPACE_DIM = 2
 
constexpr double a0 = 0.0
 
constexpr double a1 = 2.0
 
constexpr double a2 = -15.0 * a0
 
constexpr double a3 = -20.0 / 6 * a1
 
constexpr double a4 = 15.0 * a0
 
constexpr double a5 = a1
 
constexpr double a6 = -a0
 

Typedef Documentation

◆ FaceEle

Definition at line 17 of file hcurl_check_approx_in_2d.cpp.

◆ FaceEleOp

Definition at line 18 of file hcurl_check_approx_in_2d.cpp.

◆ OpDomainMass

OPerator to integrate mass matrix for least square approximation.

Examples
hcurl_check_approx_in_2d.cpp.

Definition at line 27 of file hcurl_check_approx_in_2d.cpp.

◆ OpDomainSource

using OpDomainSource = FormsIntegrators<FaceEleOp>::Assembly<PETSC>::LinearForm< GAUSS>::OpSource<BASE_DIM, BASE_DIM>

Operator to integrate the right hand side matrix for the problem.

Examples
hcurl_check_approx_in_2d.cpp.

Definition at line 34 of file hcurl_check_approx_in_2d.cpp.

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)
Examples
hcurl_check_approx_in_2d.cpp.

Definition at line 236 of file hcurl_check_approx_in_2d.cpp.

236 {
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_IMPOSIBLE_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 = smartCreateDMVector(dm);
327 auto F = smartVectorDuplicate(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}
#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_IMPOSIBLE_CASE
Definition: definitions.h:35
#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
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: DMMMoFEM.cpp:470
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
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[]
FormsIntegrators< FaceEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, BASE_DIM > OpDomainMass
OPerator to integrate mass matrix for least square approximation.
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1086
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
const double D
diffusivity
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.
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:26
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:374
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:290
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:806
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:304
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:589
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:749
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

Variable Documentation

◆ a0

constexpr double a0 = 0.0
constexpr
Examples
hcurl_check_approx_in_2d.cpp.

Definition at line 37 of file hcurl_check_approx_in_2d.cpp.

◆ a1

constexpr double a1 = 2.0
constexpr
Examples
hcurl_check_approx_in_2d.cpp.

Definition at line 38 of file hcurl_check_approx_in_2d.cpp.

◆ a2

constexpr double a2 = -15.0 * a0
constexpr
Examples
hcurl_check_approx_in_2d.cpp.

Definition at line 39 of file hcurl_check_approx_in_2d.cpp.

◆ a3

constexpr double a3 = -20.0 / 6 * a1
constexpr
Examples
hcurl_check_approx_in_2d.cpp.

Definition at line 40 of file hcurl_check_approx_in_2d.cpp.

◆ a4

constexpr double a4 = 15.0 * a0
constexpr
Examples
hcurl_check_approx_in_2d.cpp.

Definition at line 41 of file hcurl_check_approx_in_2d.cpp.

◆ a5

constexpr double a5 = a1
constexpr
Examples
hcurl_check_approx_in_2d.cpp.

Definition at line 42 of file hcurl_check_approx_in_2d.cpp.

◆ a6

constexpr double a6 = -a0
constexpr
Examples
hcurl_check_approx_in_2d.cpp.

Definition at line 43 of file hcurl_check_approx_in_2d.cpp.

◆ BASE_DIM

constexpr int BASE_DIM = 3
constexpr
Examples
hcurl_check_approx_in_2d.cpp.

Definition at line 20 of file hcurl_check_approx_in_2d.cpp.

◆ help

char help[] = "...\n\n"
static
Examples
hcurl_check_approx_in_2d.cpp.

Definition at line 15 of file hcurl_check_approx_in_2d.cpp.

◆ SPACE_DIM

constexpr int SPACE_DIM = 2
constexpr
Examples
hcurl_check_approx_in_2d.cpp.

Definition at line 21 of file hcurl_check_approx_in_2d.cpp.