v0.13.1
continuity_check_on_skeleton_with_simple_2d_for_h1.cpp
Go to the documentation of this file.
1/**
2 * \file continuity_check_on_skeleton_with_simple_2d_for_h1.cpp
3 * \ingroup mofem_simple_interface
4 \example continuity_check_on_skeleton_with_simple_2d_for_h1.cpp
5 *
6 * \brief Integration on skeleton for 2d
7 *
8 * Teting integration on skeleton and checking of continuity of hcurl space on
9 * edges.
10 */
11
12#include <MoFEM.hpp>
13using namespace MoFEM;
14
15static char help[] = "...\n\n";
16
21
22struct CommonData {
26};
27
28struct SkeletonFE : public EdgeEleOp {
29
30 struct OpFaceSide : public FaceEleOnSideOp {
31
35 elemData(elem_data) {}
36
40
41 if (type == MBEDGE && side == getEdgeSideNumber()) {
42
44
45 const double eps = 1e-12;
46 if (norm_inf(diff) > eps) {
48 LogManager::BitLineID | LogManager::BitScope);
50 MOFEM_LOG("ATOM", Sev::error)
51 << "Quad coords: " << getCoordsAtGaussPts();
52 MOFEM_LOG("ATOM", Sev::error)
53 << "Edge coords: " << getEdgeCoordsAtGaussPts();
54 MOFEM_LOG("ATOM", Sev::error) << "Difference: " << diff;
55 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
56 "Coordinates at integration pts are different");
57 }
58
59 const size_t nb_dofs = data.getN().size2();
60 const size_t nb_integration_pts = data.getN().size1();
61
62 auto t_base = data.getFTensor0N();
63 MatrixDouble *ptr_dot_elem_data = nullptr;
64 if (getFEMethod()->nInTheLoop == 0)
65 ptr_dot_elem_data = &elemData.dotEleLeft;
66 else
67 ptr_dot_elem_data = &elemData.dotEleRight;
68 MatrixDouble &dot_elem_data = *ptr_dot_elem_data;
69 dot_elem_data.resize(nb_integration_pts, nb_dofs, false);
70
71 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
72 for (size_t bb = 0; bb != nb_dofs; ++bb) {
73 dot_elem_data(gg, bb) = t_base;
74 ++t_base;
75 }
76 }
77 }
79 }
80 };
81
84
87 faceSideFe(m_field), elemData(elem_data) {
89 }
90
93
95 if (type == MBEDGE) {
96
97 const size_t nb_dofs = data.getN().size2();
98 const size_t nb_integration_pts = data.getN().size1();
99
100 MOFEM_LOG("ATOM", Sev::noisy)
101 << "Cords at integration points" << getCoordsAtGaussPts();
102
103 auto t_base = data.getFTensor0N();
104 elemData.dotEdge.resize(nb_integration_pts, nb_dofs, false);
105 elemData.dotEleLeft.resize(nb_integration_pts, 0, false);
106 elemData.dotEleRight.resize(nb_integration_pts, 0, false);
107
108 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
109 for (size_t bb = 0; bb != nb_dofs; ++bb) {
110 elemData.dotEdge(gg, bb) = t_base;
111 ++t_base;
112 }
113 }
114
116
117 auto check_continuity_of_base = [&](auto &vol_dot_data) {
119
120 if (vol_dot_data.size1() != elemData.dotEdge.size1())
121 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
122 "Inconsistent number of integration points");
123
124 if (vol_dot_data.size2() != elemData.dotEdge.size2())
125 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
126 "Inconsistent number of base functions");
127 const double eps = 1e-12;
128 for (size_t gg = 0; gg != vol_dot_data.size1(); ++gg)
129 for (size_t bb = 0; bb != vol_dot_data.size2(); ++bb) {
130 const double error =
131 std::abs(vol_dot_data(gg, bb) - elemData.dotEdge(gg, bb));
132 if (error > eps)
133 SETERRQ4(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
134 "Inconsistency (%d, %d) %3.4e != %3.4e", gg, bb,
135 vol_dot_data(gg, bb), elemData.dotEdge(gg, bb));
136 else
137 MOFEM_LOG("ATOM", Sev::noisy) << "Ok";
138 }
140 };
141
142 if (elemData.dotEleLeft.size2() != 0)
143 CHKERR check_continuity_of_base(elemData.dotEleLeft);
144 else if (elemData.dotEleRight.size2() != 0)
145 CHKERR check_continuity_of_base(elemData.dotEleRight);
146 }
148 }
149};
150
151int main(int argc, char *argv[]) {
152
153 // initialize petsc
154 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
155
156 try {
157
158 // Create MoAB database
159 moab::Core moab_core;
160 moab::Interface &moab = moab_core;
161
162 // Create MoFEM database and link it to MoAB
163 MoFEM::Core mofem_core(moab);
164 MoFEM::Interface &m_field = mofem_core;
165
166 auto core_log = logging::core::get();
167 core_log->add_sink(
168 LogManager::createSink(LogManager::getStrmSelf(), "ATOM"));
169 LogManager::setLog("ATOM");
170
171 // Register DM Manager
172 DMType dm_name = "DMMOFEM";
173 CHKERR DMRegister_MoFEM(dm_name);
174
175 // Simple interface
176 Simple *simple_interface;
177 CHKERR m_field.getInterface(simple_interface);
178 {
179 // get options from command line
180 CHKERR simple_interface->getOptions();
181 // Add ghost cells needed for evaluation on skeleton
182 simple_interface->getAddSkeletonFE() = true;
183
184 // load mesh file
185 CHKERR simple_interface->loadFile("");
186
187 auto get_base = []() -> FieldApproximationBase {
188 enum bases { AINSWORTH, DEMKOWICZ, LASTBASEOP };
189 const char *list_bases[] = {"ainsworth", "demkowicz"};
190 PetscBool flg;
191 PetscInt choice_base_value = AINSWORTH;
192 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
193 LASTBASEOP, &choice_base_value, &flg);
194 if (flg == PETSC_TRUE) {
196 if (choice_base_value == AINSWORTH)
198 else if (choice_base_value == DEMKOWICZ)
200 return base;
201 }
202 return LASTBASE;
203 };
204
205 // add fields
206 auto base = get_base();
207 CHKERR simple_interface->addDomainField("FIELD", H1, base, 1);
208 // CHKERR simple_interface->addDomainField("TEST_FIELD", L2,
209 // AINSWORTH_LEGENDRE_BASE, 1);
210 CHKERR simple_interface->addSkeletonField("FIELD", H1, base, 1);
211 // set fields order
212 CHKERR simple_interface->setFieldOrder("FIELD", 2);
213 // CHKERR simple_interface->setFieldOrder("TEST_FIELD", 1);
214 // setup problem
215 CHKERR simple_interface->setUp();
216 // get dm
217 auto dm = simple_interface->getDM();
218
219 // create elements
220 CommonData elem_data;
221 boost::shared_ptr<EdgeEle> skeleton_fe =
222 boost::shared_ptr<EdgeEle>(new EdgeEle(m_field));
223 skeleton_fe->getOpPtrVector().push_back(
224 new SkeletonFE(m_field, elem_data));
225
226 // iterate skeleton finite elements
227 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getSkeletonFEName(),
228 skeleton_fe);
229 }
230 }
232
233 // finish work cleaning memory, getting statistics, etc.
235
236 return 0;
237}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static const double eps
int main(int argc, char *argv[])
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ LASTBASE
Definition: definitions.h:69
@ 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
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ 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
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:533
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_ATTRIBUTES(channel, bit)
Add attributes to channel.
Definition: LogManager.hpp:289
#define MOFEM_LOG_FUNCTION()
Set scope.
Definition: LogManager.hpp:318
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
CoreTmp< 0 > Core
Definition: Core.hpp:1086
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
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.
MoFEMErrorCode loopSideFaces(const string fe_name, FaceElementForcesAndSourcesCoreOnSide &fe_side)
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Base face element used to integrate on skeleton.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
@ OPROW
operator doWork function is executed on FE rows
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
const std::string getSkeletonFEName() const
Get the Skeleton FE Name.
Definition: Simple.hpp:327
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 addSkeletonField(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 skeleton.
Definition: Simple.cpp:410
bool & getAddSkeletonFE()
Get the addSkeletonFE.
Definition: Simple.hpp:397
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.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
SkeletonFE(MoFEM::Interface &m_field, CommonData &elem_data)