v0.15.0
Loading...
Searching...
No Matches
tensor_divergence_operator.cpp File Reference
#include <MoFEM.hpp>

Go to the source code of this file.

Classes

struct  ElementsAndOps< DIM >
 
struct  ElementsAndOps< 2 >
 
struct  ElementsAndOps< 3 >
 

Typedefs

using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle
 
using BoundaryEle = ElementsAndOps<SPACE_DIM>::BoundaryEle
 
using AssemblyDomainEleOp = FormsIntegrators<DomainEleOp>::Assembly<A>::OpBase
 
using AssemblyBoundaryEleOp
 

Functions

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

Variables

static char help [] = "...\n\n"
 
constexpr AssemblyType A = AssemblyType::PETSC
 
constexpr IntegrationType I
 
constexpr int SPACE_DIM = EXECUTABLE_DIMENSION
 
constexpr CoordinateTypes COORD_TYPE = EXECUTABLE_COORD_TYPE
 
FTensor::Index< 'i', SPACE_DIMi
 
constexpr FieldSpace HDIV_SPACE = ElementsAndOps<SPACE_DIM>::HDIV_SPACE
 

Typedef Documentation

◆ AssemblyBoundaryEleOp

◆ AssemblyDomainEleOp

◆ BoundaryEle

◆ DomainEle

Function Documentation

◆ main()

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

Definition at line 64 of file tensor_divergence_operator.cpp.

64 {
65
66 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
67
68 try {
69
70 moab::Core mb_instance;
71 moab::Interface &moab = mb_instance;
72
73 // Create MoFEM instance
74 MoFEM::Core core(moab);
75 MoFEM::Interface &m_field = core;
76
77 // Declare elements
78 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOP };
79 const char *list_bases[] = {"ainsworth", "demkowicz"};
80 PetscBool flg;
81 PetscInt choice_base_value = AINSWORTH;
82 CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-base", list_bases,
83 LASBASETOP, &choice_base_value, &flg);
84 if (flg != PETSC_TRUE)
85 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "base not set");
87 if (choice_base_value == AINSWORTH)
89 else if (choice_base_value == DEMKOWICZ)
91 int order = 4;
92 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order,
93 PETSC_NULLPTR);
94
95 // Register DM Manager
96 DMType dm_name = "DMMOFEM";
97 CHKERR DMRegister_MoFEM(dm_name);
98
99 // Add logging channel for example
100 auto core_log = logging::core::get();
101 core_log->add_sink(
103 LogManager::setLog("ATOM");
104 MOFEM_LOG_TAG("ATOM", "ATOM");
105
106 // Simple interface
107 auto simple = m_field.getInterface<Simple>();
108
109 // get options from command line
111 // load mesh file
112 CHKERR simple->loadFile();
113
114 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
115 CHKERR simple->addDomainField("SIGMA", HDIV_SPACE, DEMKOWICZ_JACOBI_BASE,
116 SPACE_DIM);
117 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
118 CHKERR simple->addBoundaryField("SIGMA", HDIV_SPACE, DEMKOWICZ_JACOBI_BASE,
119 SPACE_DIM);
120
121 CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
122
123 // set fields order, i.e. for most first cases order is sufficient.
124 CHKERR simple->setFieldOrder("U", order);
125 // CHKERR simple->setFieldOrder("SIGMA", order);
126 CHKERR simple->setFieldOrder("GEOMETRY", 2);
127
128 auto get_skin = [&]() {
129 Range body_ents;
130 CHKERR m_field.get_moab().get_entities_by_dimension(0, SPACE_DIM,
131 body_ents);
132 Skinner skin(&m_field.get_moab());
133 Range skin_ents;
134 CHKERR skin.find_skin(0, body_ents, false, skin_ents);
135 return skin_ents;
136 };
137
138 auto filter_true_skin = [&](auto skin) {
139 Range boundary_ents;
140 ParallelComm *pcomm =
141 ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
142 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
143 PSTATUS_NOT, -1, &boundary_ents);
144 return boundary_ents;
145 };
146
147 auto boundary_ents = filter_true_skin(get_skin());
148 CHKERR simple->setFieldOrder("SIGMA", 0);
149 CHKERR simple->setFieldOrder("SIGMA", order, &boundary_ents);
150
151 // setup problem
152 CHKERR simple->setUp();
153
154 auto bc_mng = m_field.getInterface<BcManager>();
155 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM",
156 "U", 0, MAX_DOFS_ON_ENTITY, true);
157 CHKERR bc_mng->removeBlockDOFsOnEntities(
158 simple->getProblemName(), "SYM", "SIGMA", 0, MAX_DOFS_ON_ENTITY, true);
159
160 auto project_ho_geometry = [&]() {
161 Projection10NodeCoordsOnField ent_method(m_field, "GEOMETRY");
162 return m_field.loop_dofs("GEOMETRY", ent_method);
163 };
164 CHKERR project_ho_geometry();
165
166 // get operators tester
167 auto opt = m_field.getInterface<OperatorsTester>(); // get interface to
168 // OperatorsTester
169 auto pip_mng = m_field.getInterface<PipelineManager>(); // get interface to
170 // pipeline manager
171
172 pip_mng->getOpDomainLhsPipeline().clear();
173 pip_mng->getOpDomainRhsPipeline().clear();
174
175 // integration rule
176 auto rule = [&](int, int, int p) { return 2 * p + 2; };
177 CHKERR pip_mng->setDomainRhsIntegrationRule(rule);
178 CHKERR pip_mng->setDomainLhsIntegrationRule(rule);
179 CHKERR pip_mng->setBoundaryRhsIntegrationRule(rule);
180
181 auto x = opt->setRandomFields(simple->getDM(),
182 {{"U", {-1, 1}}, {"SIGMA", {-1, 1}}});
183 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), x, INSERT_VALUES,
184 SCATTER_REVERSE);
185
186 // set integration rules
187
188 auto jacobian = [&](double r) {
189 if constexpr (COORD_TYPE == CYLINDRICAL)
190 return 2 * M_PI * r;
191 else
192 return 1.;
193 };
194
195 auto beta_domain = [jacobian](double r, double, double) {
196 return jacobian(r);
197 };
198 auto beta_bdy = [jacobian](double r, double, double) {
199 return -jacobian(r);
200 };
201
202 auto post_proc = [&](auto dm, auto f_res, auto out_name) {
204 auto post_proc_fe =
205 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(m_field);
206
208
209 auto sigma_ptr = boost::make_shared<MatrixDouble>();
210 post_proc_fe->getOpPtrVector().push_back(
212 sigma_ptr));
213 auto u_ptr = boost::make_shared<MatrixDouble>();
214 post_proc_fe->getOpPtrVector().push_back(
216
217 post_proc_fe->getOpPtrVector().push_back(
218
219 new OpPPMap(
220
221 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
222 {}, {{"U", u_ptr}}, {{"SIGMA", sigma_ptr}}, {})
223
224 );
225
226 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), f_res, INSERT_VALUES,
227 SCATTER_REVERSE);
228 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
229 post_proc_fe);
230 post_proc_fe->writeFile(out_name);
232 };
233
234 auto test_consistency_of_domain_and_bdy_integrals = [&]() {
236 using OpMixDivURhs =
238 GAUSS>::OpMixDivTimesU<3, SPACE_DIM, SPACE_DIM, COORD_TYPE>;
239 using OpMixLambdaGradURhs = FormsIntegrators<DomainEleOp>::Assembly<
241 using OpMixLambdaGradURhs = FormsIntegrators<DomainEleOp>::Assembly<
243 using OpMixUTimesDivLambdaRhs = FormsIntegrators<DomainEleOp>::Assembly<
245
246 using OpMixUTimesLambdaRhs = FormsIntegrators<DomainEleOp>::Assembly<
248
249 using OpMixNormalLambdaURhs = FormsIntegrators<BoundaryEleOp>::Assembly<
251 using OpUTimeTractionRhs = FormsIntegrators<BoundaryEleOp>::Assembly<
253
254 auto ops_rhs_interior = [&](auto &pip) {
257 "GEOMETRY");
258 auto u_ptr = boost::make_shared<MatrixDouble>();
259 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
260 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
262 "U", grad_u_ptr));
263
264 pip.push_back(new OpMixDivURhs("SIGMA", u_ptr, beta_domain));
265 pip.push_back(
266 new OpMixLambdaGradURhs("SIGMA", grad_u_ptr, beta_domain));
267
268 auto sigma_ptr = boost::make_shared<MatrixDouble>();
269 auto sigma_div_ptr = boost::make_shared<MatrixDouble>();
271 COORD_TYPE>(
272 "SIGMA", sigma_div_ptr));
274 "SIGMA", sigma_ptr));
275 pip.push_back(
276 new OpMixUTimesDivLambdaRhs("U", sigma_div_ptr, beta_domain));
277 pip.push_back(new OpMixUTimesLambdaRhs("U", sigma_ptr, beta_domain));
278
280 };
281
282 auto ops_rhs_boundary = [&](auto &pip) {
285 "GEOMETRY");
286 auto u_ptr = boost::make_shared<MatrixDouble>();
287 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
288 auto traction_ptr = boost::make_shared<MatrixDouble>();
290 "SIGMA", traction_ptr));
291
292 // We have to integrate on curved face geometry, thus integration weight
293 // have to adjusted.
294 pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
295 pip.push_back(new OpMixNormalLambdaURhs("SIGMA", u_ptr, beta_bdy));
296 pip.push_back(new OpUTimeTractionRhs("U", traction_ptr, beta_bdy));
298 };
299
300 CHKERR ops_rhs_interior(pip_mng->getOpDomainRhsPipeline());
301 CHKERR ops_rhs_boundary(pip_mng->getOpBoundaryRhsPipeline());
302
303 auto f = createDMVector(simple->getDM());
304 pip_mng->getDomainRhsFE()->f = f;
305 pip_mng->getBoundaryRhsFE()->f = f;
306
307 CHKERR VecZeroEntries(f);
308
310 simple->getDomainFEName(),
311 pip_mng->getDomainRhsFE());
313 simple->getBoundaryFEName(),
314 pip_mng->getBoundaryRhsFE());
315
316 CHKERR VecAssemblyBegin(f);
317 CHKERR VecAssemblyEnd(f);
318
319 double f_nrm2;
320 CHKERR VecNorm(f, NORM_2, &f_nrm2);
321
322 MOFEM_LOG("ATOM", Sev::inform) << "f_norm2 = " << f_nrm2;
323 if (std::fabs(f_nrm2) > 1e-10) {
324 CHKERR post_proc(simple->getDM(), f,
325 "tensor_divergence_operator_res_vec.h5m");
326 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "Test failed");
327 }
328
330 };
331
332 // Testing Lhs domain operators
333
334 auto test_lhs_ops = [&]() {
336 using OpMixDivULhs = FormsIntegrators<DomainEleOp>::Assembly<
338 using OpLambdaGraULhs = FormsIntegrators<DomainEleOp>::Assembly<
340
341 auto op_lhs_domain = [&](auto &pip) {
344 "GEOMETRY");
345
346 auto unity = []() { return 1; };
347 pip.push_back(
348 new OpMixDivULhs("SIGMA", "U", unity, beta_domain, true, false));
349 pip.push_back(
350 new OpLambdaGraULhs("SIGMA", "U", unity, beta_domain, true, false));
352 };
353
354 CHKERR op_lhs_domain(pip_mng->getOpDomainLhsPipeline());
355
356 auto diff_x = opt->setRandomFields(simple->getDM(),
357 {{"U", {-1, 1}}, {"SIGMA", {-1, 1}}});
358 constexpr double eps = 1e-5;
359 auto diff_res = opt->checkCentralFiniteDifference(
360 simple->getDM(), simple->getDomainFEName(), pip_mng->getDomainRhsFE(),
361 pip_mng->getDomainLhsFE(), x, SmartPetscObj<Vec>(),
362 SmartPetscObj<Vec>(), diff_x, 0, 1, eps);
363
364 // Calculate norm of difference between directive calculated from finite
365 // difference, and tangent matrix.
366 double fnorm_res;
367 CHKERR VecNorm(diff_res, NORM_2, &fnorm_res);
368 MOFEM_LOG_C("ATOM", Sev::inform, "Test Lhs OPs %3.4e", fnorm_res);
369 if (std::fabs(fnorm_res) > 1e-8)
370 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "Test failed");
371
373 };
374
375 CHKERR test_consistency_of_domain_and_bdy_integrals();
376 CHKERR test_lhs_ops();
377 }
379
381}
static auto filter_true_skin(MoFEM::Interface &m_field, Range &&skin)
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
static const double eps
#define CATCH_ERRORS
Catch errors.
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
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
@ H1
continuous field
Definition definitions.h:85
@ HDIV
field with continuous normal traction
Definition definitions.h:87
#define MYPCOMM_INDEX
default communicator number PCOMM
#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
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
@ CYLINDRICAL
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
constexpr int order
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:514
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
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:576
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1102
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
int r
Definition sdf.py:8
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition BcManager.hpp:29
MoFEMErrorCode removeBlockDOFsOnEntities(const std::string problem_name, const std::string block_name, const std::string field_name, int lo, int hi, bool get_low_dim_ents=true, bool is_distributed_mesh=true)
Remove DOFs from problem.
Definition BcManager.cpp:72
virtual moab::Interface & get_moab()=0
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:118
Deprecated interface functions.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Calculate divergence of tonsorial field using vectorial base.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Calculate trace of vector (Hdiv/Hcurl) space.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Get values at integration pts for tensor field rank 1, i.e. vector field.
Post post-proc data at points from hash maps.
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
SmartPetscObj< Vec > setRandomFields(SmartPetscObj< DM > dm, std::vector< RandomFieldData > random_fields, boost::shared_ptr< Range > ents=nullptr)
Generate random fields.
PipelineManager interface.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
static char help[]
constexpr int SPACE_DIM
constexpr FieldSpace HDIV_SPACE
constexpr CoordinateTypes COORD_TYPE

Variable Documentation

◆ A

AssemblyType A = AssemblyType::PETSC
constexpr

Definition at line 31 of file tensor_divergence_operator.cpp.

◆ COORD_TYPE

CoordinateTypes COORD_TYPE = EXECUTABLE_COORD_TYPE
constexpr

◆ HDIV_SPACE

FieldSpace HDIV_SPACE = ElementsAndOps<SPACE_DIM>::HDIV_SPACE
constexpr

◆ help

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

Definition at line 29 of file tensor_divergence_operator.cpp.

◆ I

IntegrationType I
constexpr
Initial value:
=
IntegrationType::GAUSS

Definition at line 32 of file tensor_divergence_operator.cpp.

◆ i

Definition at line 38 of file tensor_divergence_operator.cpp.

◆ SPACE_DIM

int SPACE_DIM = EXECUTABLE_DIMENSION
constexpr

Definition at line 35 of file tensor_divergence_operator.cpp.