v0.14.0
Loading...
Searching...
No Matches
testing_jacobian_of_surface_pressure_element.cpp File Reference

Go to the source code of this file.

Functions

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

Variables

static char help [] = "\n"
 

Function Documentation

◆ main()

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

Definition at line 19 of file testing_jacobian_of_surface_pressure_element.cpp.

19 {
20
21 // Initialize MoFEM
22 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
23
24 // Create mesh database
25 moab::Core mb_instance; // create database
26 moab::Interface &moab = mb_instance; // create interface to database
27
28 try {
29 // Create MoFEM database and link it to MoAB
30 MoFEM::Core core(moab);
31 MoFEM::Interface &m_field = core;
32
33 PetscInt order_x = 2;
34 PetscInt order_X = 2;
35 PetscBool flg = PETSC_TRUE;
36
37 PetscBool test_jacobian = PETSC_FALSE;
38 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test_jacobian", &test_jacobian,
39 PETSC_NULL);
40 CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-order_spat", &order_x,
41 &flg);
42 CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-order_mat", &order_X,
43 &flg);
44
45 CHKERR DMRegister_MoFEM("DMMOFEM");
46
47 Simple *si = m_field.getInterface<MoFEM::Simple>();
48
49 CHKERR si->getOptions();
50 CHKERR si->loadFile();
51 CHKERR si->addDomainField("SPATIAL_POSITION", H1, AINSWORTH_LOBATTO_BASE,
52 3);
53 CHKERR si->addBoundaryField("SPATIAL_POSITION", H1, AINSWORTH_LOBATTO_BASE,
54 3);
55 CHKERR si->setFieldOrder("SPATIAL_POSITION", order_x);
56
57 CHKERR si->addDomainField("MESH_NODE_POSITIONS", H1,
59 CHKERR si->addBoundaryField("MESH_NODE_POSITIONS", H1,
61 CHKERR si->setFieldOrder("MESH_NODE_POSITIONS", order_X);
62
63 Range triangle_springs;
65 if (it->getName().compare(0, 9, "SPRING_BC") == 0) {
66 CHKERR m_field.get_moab().get_entities_by_type(it->meshset, MBTRI,
67 triangle_springs, true);
68 }
69 }
70
71 // Add spring boundary condition applied on surfaces.
72 CHKERR MetaSpringBC::addSpringElements(m_field, "SPATIAL_POSITION",
73 "MESH_NODE_POSITIONS");
75 m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS", triangle_springs);
76 si->getOtherFiniteElements().push_back("SPRING");
77 si->getOtherFiniteElements().push_back("SPRING_ALE");
78 CHKERR si->setUp();
79
80 // create DM
81 DM dm;
82 CHKERR si->getDM(&dm);
83
84 PetscRandom rctx;
85 PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
86
87 auto set_coord = [&](VectorAdaptor &&field_data, double *x, double *y,
88 double *z) {
90 double value;
91 double scale = 0.5;
92 PetscRandomGetValue(rctx, &value);
93 field_data[0] = (*x) + (value - 0.5) * scale;
94 PetscRandomGetValue(rctx, &value);
95 field_data[1] = (*y) + (value - 0.5) * scale;
96 PetscRandomGetValue(rctx, &value);
97 field_data[2] = (*z) + (value - 0.5) * scale;
99 };
100
101 CHKERR m_field.getInterface<FieldBlas>()->setVertexDofs(set_coord,
102 "SPATIAL_POSITION");
103 CHKERR m_field.getInterface<FieldBlas>()->setVertexDofs(
104 set_coord, "MESH_NODE_POSITIONS");
105
106 PetscRandomDestroy(&rctx);
107
108 boost::shared_ptr<NeumannForcesSurface> surfacePressure(
109 new NeumannForcesSurface(m_field));
110
111 boost::shared_ptr<NeumannForcesSurface::MyTriangleFE> fe_rhs_ptr(
112 surfacePressure, &(surfacePressure->getLoopFe()));
113 boost::shared_ptr<NeumannForcesSurface::MyTriangleFE> fe_lhs_ptr(
114 surfacePressure, &(surfacePressure->getLoopFeLhs()));
115 boost::shared_ptr<NeumannForcesSurface::MyTriangleFE> fe_mat_rhs_ptr(
116 surfacePressure, &(surfacePressure->getLoopFeMatRhs()));
117 boost::shared_ptr<NeumannForcesSurface::MyTriangleFE> fe_mat_lhs_ptr(
118 surfacePressure, &(surfacePressure->getLoopFeMatLhs()));
119
120 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", *fe_rhs_ptr, false, false);
121 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", *fe_lhs_ptr, false, false);
122 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", *fe_mat_rhs_ptr, false, false);
123 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", *fe_mat_lhs_ptr, false, false);
124
125 boost::shared_ptr<NeumannForcesSurface> surfaceForce(
126 new NeumannForcesSurface(m_field));
127
128 boost::shared_ptr<NeumannForcesSurface::MyTriangleFE>
129 fe_rhs_surface_force_ptr(surfaceForce, &(surfaceForce->getLoopFe()));
130 boost::shared_ptr<NeumannForcesSurface::MyTriangleFE>
131 fe_lhs_surface_force_ptr(surfaceForce, &(surfaceForce->getLoopFeLhs()));
132 boost::shared_ptr<NeumannForcesSurface::MyTriangleFE>
133 fe_mat_rhs_surface_force_ptr(surfaceForce,
134 &(surfaceForce->getLoopFeMatRhs()));
135 boost::shared_ptr<NeumannForcesSurface::MyTriangleFE>
136 fe_mat_lhs_surface_force_ptr(surfaceForce,
137 &(surfaceForce->getLoopFeMatLhs()));
138
139 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", *fe_rhs_surface_force_ptr,
140 false, false);
141 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", *fe_lhs_surface_force_ptr,
142 false, false);
143 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", *fe_mat_rhs_surface_force_ptr,
144 false, false);
145 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", *fe_mat_lhs_surface_force_ptr,
146 false, false);
147
148 Range nodes;
149 CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes, false);
150
151 nodes.pop_front();
152 nodes.pop_back();
153
154 boost::shared_ptr<NeumannForcesSurface::DataAtIntegrationPts> dataAtPts =
155 boost::make_shared<NeumannForcesSurface::DataAtIntegrationPts>();
156
157 dataAtPts->forcesOnlyOnEntitiesRow = nodes;
158
160 if (bit->getName().compare(0, 8, "PRESSURE") == 0) {
161 CHKERR surfacePressure->addPressure("SPATIAL_POSITION", PETSC_NULL,
162 bit->getMeshsetId(), true, true);
163 CHKERR surfacePressure->addPressureAle(
164 "SPATIAL_POSITION", "MESH_NODE_POSITIONS", dataAtPts,
165 si->getDomainFEName(), PETSC_NULL, PETSC_NULL, bit->getMeshsetId(),
166 true, true);
167 }
168 }
169
170 const string block_set_force_name("FORCE");
171 // search for block named FORCE and add its attributes to FORCE_FE element
173 bit)) {
174 CHKERR surfaceForce->addForce("SPATIAL_POSITION", PETSC_NULL,
175 (bit->getMeshsetId()), true, false);
176 CHKERR surfaceForce->addForceAle(
177 "SPATIAL_POSITION", "MESH_NODE_POSITIONS", dataAtPts,
178 si->getDomainFEName(), PETSC_NULL, PETSC_NULL, bit->getMeshsetId(),
179 true, false, false);
180 }
181
183 if (it->getName().compare(0, block_set_force_name.length(),
184 block_set_force_name) == 0) {
185 CHKERR surfaceForce->addForce("SPATIAL_POSITION", PETSC_NULL,
186 (it->getMeshsetId()), true, true);
187 CHKERR surfaceForce->addForceAle(
188 "SPATIAL_POSITION", "MESH_NODE_POSITIONS", dataAtPts,
189 si->getDomainFEName(), PETSC_NULL, PETSC_NULL, it->getMeshsetId(),
190 true, true, false);
191 }
192 }
193
194 // Implementation of spring element
195 // Create new instances of face elements for springs
196 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
198 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
201 m_field, fe_spring_lhs_ptr, fe_spring_rhs_ptr, "SPATIAL_POSITION",
202 "MESH_NODE_POSITIONS");
203
204 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ale_ptr_dx(
206 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ale_ptr_dX(
208
209 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ale_ptr(
211
212 boost::shared_ptr<MetaSpringBC::DataAtIntegrationPtsSprings>
213 data_at_spring_gp =
214 boost::make_shared<MetaSpringBC::DataAtIntegrationPtsSprings>(
215 m_field);
216
217 Range spring_ale_nodes;
218 CHKERR moab.get_connectivity(triangle_springs, spring_ale_nodes, true);
219
220 data_at_spring_gp->forcesOnlyOnEntitiesRow = spring_ale_nodes;
221
223 m_field, fe_spring_lhs_ale_ptr_dx, fe_spring_lhs_ale_ptr_dX,
224 fe_spring_rhs_ale_ptr, data_at_spring_gp, "SPATIAL_POSITION",
225 "MESH_NODE_POSITIONS", si->getDomainFEName());
226
227 CHKERR DMMoFEMSNESSetJacobian(dm, "SPRING", fe_spring_lhs_ptr, PETSC_NULL,
228 PETSC_NULL);
229 CHKERR DMMoFEMSNESSetFunction(dm, "SPRING", fe_spring_rhs_ptr, PETSC_NULL,
230 PETSC_NULL);
231
232 CHKERR DMMoFEMSNESSetJacobian(dm, "SPRING_ALE", fe_spring_lhs_ale_ptr_dx,
233 PETSC_NULL, PETSC_NULL);
234 CHKERR DMMoFEMSNESSetJacobian(dm, "SPRING_ALE", fe_spring_lhs_ale_ptr_dX,
235 PETSC_NULL, PETSC_NULL);
236 CHKERR DMMoFEMSNESSetFunction(dm, "SPRING_ALE", fe_spring_rhs_ale_ptr,
237 PETSC_NULL, PETSC_NULL);
238
239 CHKERR DMMoFEMSNESSetJacobian(dm, si->getBoundaryFEName(), fe_lhs_ptr,
240 nullptr, nullptr);
241 CHKERR DMMoFEMSNESSetFunction(dm, si->getBoundaryFEName(), fe_rhs_ptr,
242 nullptr, nullptr);
243
244 CHKERR DMMoFEMSNESSetJacobian(dm, si->getBoundaryFEName(), fe_mat_lhs_ptr,
245 nullptr, nullptr);
246 CHKERR DMMoFEMSNESSetFunction(dm, si->getBoundaryFEName(), fe_mat_rhs_ptr,
247 nullptr, nullptr);
248
249 // Surface force ALE
251 fe_lhs_surface_force_ptr, nullptr, nullptr);
253 fe_rhs_surface_force_ptr, nullptr, nullptr);
254
256 fe_mat_lhs_surface_force_ptr, nullptr,
257 nullptr);
259 fe_mat_rhs_surface_force_ptr, nullptr,
260 nullptr);
261
262 Vec x, f;
263 CHKERR DMCreateGlobalVector(dm, &x);
264 CHKERR VecDuplicate(x, &f);
265 CHKERR VecSetOption(f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
266 CHKERR DMoFEMMeshToLocalVector(dm, x, INSERT_VALUES, SCATTER_FORWARD);
267
268 Mat A, fdA;
269 CHKERR DMCreateMatrix(dm, &A);
270 CHKERR MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &fdA);
271
272 if (test_jacobian == PETSC_TRUE) {
273 char testing_options[] =
274 "-snes_test_jacobian -snes_test_jacobian_display "
275 "-snes_no_convergence_test -snes_atol 0 -snes_rtol 0 -snes_max_it 1 ";
276 //"-pc_type none";
277 CHKERR PetscOptionsInsertString(NULL, testing_options);
278 } else {
279 char testing_options[] = "-snes_no_convergence_test -snes_atol 0 "
280 "-snes_rtol 0 "
281 "-snes_max_it 1 ";
282 //"-pc_type none";
283 CHKERR PetscOptionsInsertString(NULL, testing_options);
284 }
285
286 SNES snes;
287 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
288 MoFEM::SnesCtx *snes_ctx;
289 CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
290 CHKERR SNESSetFunction(snes, f, SnesRhs, snes_ctx);
291 CHKERR SNESSetJacobian(snes, A, A, SnesMat, snes_ctx);
292 CHKERR SNESSetFromOptions(snes);
293
294 CHKERR SNESSolve(snes, NULL, x);
295
296 if (test_jacobian == PETSC_FALSE) {
297 double nrm_A0;
298 CHKERR MatNorm(A, NORM_INFINITY, &nrm_A0);
299
300 char testing_options_fd[] = "-snes_fd";
301 CHKERR PetscOptionsInsertString(NULL, testing_options_fd);
302
303 CHKERR SNESSetFunction(snes, f, SnesRhs, snes_ctx);
304 CHKERR SNESSetJacobian(snes, fdA, fdA, SnesMat, snes_ctx);
305 CHKERR SNESSetFromOptions(snes);
306
307 CHKERR SNESSolve(snes, NULL, x);
308 CHKERR MatAXPY(A, -1, fdA, SUBSET_NONZERO_PATTERN);
309
310 double nrm_A;
311 CHKERR MatNorm(A, NORM_INFINITY, &nrm_A);
312 PetscPrintf(PETSC_COMM_WORLD, "Matrix norms %3.4e %3.4e\n", nrm_A,
313 nrm_A / nrm_A0);
314 nrm_A /= nrm_A0;
315
316 const double tol = 1e-7;
317 if (nrm_A > tol) {
318 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
319 "Difference between hand-calculated tangent matrix and finite "
320 "difference matrix is too big");
321 }
322 }
323
324 CHKERR VecDestroy(&x);
325 CHKERR VecDestroy(&f);
326 CHKERR MatDestroy(&A);
327 CHKERR MatDestroy(&fdA);
328 CHKERR SNESDestroy(&snes);
329
330 // destroy DM
331 CHKERR DMDestroy(&dm);
332 }
334
335 // finish work cleaning memory, getting statistics, etc
337
338 return 0;
339}
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
@ AINSWORTH_LOBATTO_BASE
Definition definitions.h:62
@ 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 ...
@ FORCESET
@ NODESET
@ BLOCKSET
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
Definition DMMoFEM.cpp:718
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 DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
Definition DMMoFEM.cpp:759
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition DMMoFEM.cpp:1094
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
auto bit
set bit
double tol
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition Types.hpp:115
PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
Definition SnesCtx.cpp:139
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
Definition SnesCtx.cpp:27
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
constexpr AssemblyType A
double scale
Definition plastic.cpp:119
static MoFEMErrorCode setSpringOperators(MoFEM::Interface &m_field, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_lhs_ptr, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_rhs_ptr, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", double stiffness_scale=1.)
Implementation of spring element. Set operators to calculate LHS and RHS.
static MoFEMErrorCode addSpringElementsALE(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions, Range &spring_triangles)
Declare spring element in ALE.
static MoFEMErrorCode addSpringElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare spring element.
static MoFEMErrorCode setSpringOperatorsMaterial(MoFEM::Interface &m_field, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_lhs_ptr_dx, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_lhs_ptr_dX, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_rhs_ptr, boost::shared_ptr< DataAtIntegrationPtsSprings > data_at_integration_pts, const std::string field_name, const std::string mesh_nodals_positions, std::string side_fe_name)
Implementation of spring element. Set operators to calculate LHS and RHS.
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:112
Deprecated interface functions.
Basic algebra on fields.
Definition FieldBlas.hpp:21
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
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
Definition Simple.hpp:394
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition Simple.cpp:194
MoFEMErrorCode addBoundaryField(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 boundary.
Definition Simple.cpp:358
std::vector< std::string > & getOtherFiniteElements()
Get the Other Finite Elements.
Definition Simple.hpp:462
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
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition Simple.hpp:387
Interface for nonlinear (SNES) solver.
Definition SnesCtx.hpp:13
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Finite element and operators to apply force/pressures applied to surfaces.

Variable Documentation

◆ help

char help[] = "\n"
static