v0.14.0
Loading...
Searching...
No Matches
continuity_check_on_skeleton_with_simple_2d_for_hcurl.cpp
Go to the documentation of this file.
1/**
2 * \file continuity_check_on_skeleton_with_simple_2d_for_hcurl.cpp
3 * \ingroup mofem_simple_interface
4 * \example continuity_check_on_skeleton_with_simple_2d_for_hcurl.cpp
5 *
6 * \brief Integration on skeleton for 2d
7 *
8 * Testing 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
34 : FaceEleOnSideOp("FIELD", UserDataOperator::OPROW), elemData(elem_data) {}
35
39
40 if (type == MBEDGE && side == getEdgeSideNumber()) {
41
43
44 const double eps = 1e-12;
45 if (norm_inf(diff) > eps)
46 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
47 "coordinates at integration pts are different");
48
49 const size_t nb_dofs = data.getN().size2() / 3;
50 const size_t nb_integration_pts = data.getN().size1();
51
52 auto t_tangent = getFTensor1Direction();
53 auto t_hcurl_base = data.getFTensor1N<3>();
55 MatrixDouble *ptr_dot_elem_data = nullptr;
56 if (getFEMethod()->nInTheLoop == 0)
57 ptr_dot_elem_data = &elemData.dotEleLeft;
58 else
59 ptr_dot_elem_data = &elemData.dotEleRight;
60 MatrixDouble &dot_elem_data = *ptr_dot_elem_data;
61 dot_elem_data.resize(nb_integration_pts, nb_dofs, false);
62
63 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
64 for (size_t bb = 0; bb != nb_dofs; ++bb) {
65 dot_elem_data(gg, bb) = t_tangent(i) * t_hcurl_base(i);
66 ++t_hcurl_base;
67 }
68 }
69 }
71 }
72 };
73
76
79 faceSideFe(m_field), elemData(elem_data) {
82 }
83
86
88 if (type == MBEDGE) {
89
90 const size_t nb_dofs = data.getN().size2() / 3;
91 const size_t nb_integration_pts = data.getN().size1();
92
93 auto t_tangent = getFTensor1Direction();
94 auto t_hcurl_base = data.getFTensor1N<3>();
96 elemData.dotEdge.resize(nb_integration_pts, nb_dofs, false);
97 elemData.dotEleLeft.resize(nb_integration_pts, 0, false);
98 elemData.dotEleRight.resize(nb_integration_pts, 0, false);
99
100 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
101 for (size_t bb = 0; bb != nb_dofs; ++bb) {
102 elemData.dotEdge(gg, bb) = t_tangent(i) * t_hcurl_base(i);
103 ++t_hcurl_base;
104 }
105 }
106
108
109 auto check_continuity_of_base = [&](auto &vol_dot_data) {
111
112 if (vol_dot_data.size1() != elemData.dotEdge.size1())
113 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
114 "Inconsistent number of integration points");
115
116 if (vol_dot_data.size2() != elemData.dotEdge.size2())
117 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
118 "Inconsistent number of base functions");
119 const double eps = 1e-12;
120 for (size_t gg = 0; gg != vol_dot_data.size1(); ++gg)
121 for (size_t bb = 0; bb != vol_dot_data.size2(); ++bb) {
122 const double error =
123 std::abs(vol_dot_data(gg, bb) - elemData.dotEdge(gg, bb));
124 if (error > eps)
125 SETERRQ4(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
126 "Inconsistency (%d, %d) %3.4e != %3.4e", gg, bb,
127 vol_dot_data(gg, bb), elemData.dotEdge(gg, bb));
128 else
129 MOFEM_LOG("ATOM", Sev::noisy) << "Ok";
130
131 }
133 };
134
135 if (elemData.dotEleLeft.size2() != 0)
136 CHKERR check_continuity_of_base(elemData.dotEleLeft);
137 else if (elemData.dotEleRight.size2() != 0)
138 CHKERR check_continuity_of_base(elemData.dotEleRight);
139
140 }
142 }
143};
144
145int main(int argc, char *argv[]) {
146
147 // initialize petsc
148 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
149
150 try {
151
152 // Create MoAB database
153 moab::Core moab_core;
154 moab::Interface &moab = moab_core;
155
156 // Create MoFEM database and link it to MoAB
157 MoFEM::Core mofem_core(moab);
158 MoFEM::Interface &m_field = mofem_core;
159
160 auto core_log = logging::core::get();
161 core_log->add_sink(
163 LogManager::setLog("ATOM");
164
165 // Register DM Manager
166 DMType dm_name = "DMMOFEM";
167 CHKERR DMRegister_MoFEM(dm_name);
168
169 // Simple interface
170 Simple *simple_interface;
171 CHKERR m_field.getInterface(simple_interface);
172 {
173 // get options from command line
174 CHKERR simple_interface->getOptions();
175 // load mesh file
176 CHKERR simple_interface->loadFile("");
177
178 auto get_base = []() -> FieldApproximationBase {
179 enum bases { AINSWORTH, DEMKOWICZ, LASTBASEOP };
180 const char *list_bases[] = {"ainsworth", "demkowicz"};
181 PetscBool flg;
182 PetscInt choice_base_value = AINSWORTH;
183 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
184 LASTBASEOP, &choice_base_value, &flg);
185 if (flg == PETSC_TRUE) {
187 if (choice_base_value == AINSWORTH)
189 else if (choice_base_value == DEMKOWICZ)
191 return base;
192 }
193 return LASTBASE;
194 };
195
196 // add fields
197 auto base = get_base();
198 CHKERR simple_interface->addDomainField("FIELD", HCURL, base, 1);
199 CHKERR simple_interface->addSkeletonField("FIELD", HCURL, base, 1);
200 // set fields order
201 CHKERR simple_interface->setFieldOrder("FIELD", 3);
202 // setup problem
203 CHKERR simple_interface->setUp();
204 // get dm
205 auto dm = simple_interface->getDM();
206
207 // create elements
208 CommonData elem_data;
209 boost::shared_ptr<EdgeEle> skeleton_fe =
210 boost::shared_ptr<EdgeEle>(new EdgeEle(m_field));
211
212 CHKERR AddHOOps<1, 2, 2>::add(skeleton_fe->getOpPtrVector(), {HCURL});
213 skeleton_fe->getOpPtrVector().push_back(
214 new SkeletonFE(m_field, elem_data));
215
216 // iterate skeleton finite elements
217 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getSkeletonFEName(),
218 skeleton_fe);
219 }
220 }
222
223 // finish work cleaning memory, getting statistics, etc.
225
226 return 0;
227}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
int main()
Definition: adol-c_atom.cpp:46
static const double eps
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
@ 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_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: DMMoFEM.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: DMMoFEM.cpp:572
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
FTensor::Index< 'i', SPACE_DIM > i
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Add operators pushing bases from local to physical configuration.
int nInTheLoop
number currently of processed method
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)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
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_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
static boost::shared_ptr< std::ostream > getStrmSelf()
Get the strm self object.
Definition: LogManager.cpp:340
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
const std::string getSkeletonFEName() const
Get the Skeleton FE Name.
Definition: Simple.hpp:355
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 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:300
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.
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)