v0.14.0
Loading...
Searching...
No Matches
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 OpDomainMass
 OPerator to integrate mass matrix for least square approximation.
 
using OpDomainSource
 Operator to integrate the right hand side matrix for the problem.
 

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

◆ OpDomainMass

using OpDomainMass
Initial value:
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition seepage.cpp:56

OPerator to integrate mass matrix for least square approximation.

Definition at line 27 of file hcurl_check_approx_in_2d.cpp.

◆ OpDomainSource

Initial value:

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

Definition at line 34 of file hcurl_check_approx_in_2d.cpp.

Function Documentation

◆ main()

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

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_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}
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
Definition definitions.h:58
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
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 ...
@ MOFEM_IMPOSSIBLE_CASE
Definition definitions.h:35
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
constexpr int order
@ 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:523
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1099
FormsIntegrators< FaceEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, BASE_DIM > OpDomainMass
OPerator to integrate mass matrix for least square approximation.
static char help[]
FormsIntegrators< FaceEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, BASE_DIM > OpDomainSource
Operator to integrate the right hand side matrix for the problem.
double D
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
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:27
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 loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition Simple.cpp:194
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition Simple.cpp:829
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition Simple.cpp:599
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition Simple.cpp:765
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.

Variable Documentation

◆ a0

◆ a1

◆ a2

◆ a3

double a3 = -20.0 / 6 * a1
constexpr

◆ a4

◆ a5

double a5 = a1
constexpr

◆ a6

double a6 = -a0
constexpr

◆ BASE_DIM

int BASE_DIM = 3
constexpr

Definition at line 20 of file hcurl_check_approx_in_2d.cpp.

◆ help

char help[] = "...\n\n"
static

Definition at line 15 of file hcurl_check_approx_in_2d.cpp.

◆ SPACE_DIM

int SPACE_DIM = 2
constexpr

Definition at line 21 of file hcurl_check_approx_in_2d.cpp.