v0.16.0
Loading...
Searching...
No Matches
MagneticElement.hpp
Go to the documentation of this file.
1/** \file MagneticElement.hpp
2 * \brief Implementation of magnetic element
3 * \ingroup maxwell_element
4 * \example mofem/tutorials/max-0_magnetostatics/src/MagneticElement.hpp
5 *
6 */
7
8#ifndef __MAGNETICELEMENT_HPP__
9#define __MAGNETICELEMENT_HPP__
10
11/**
12 * \brief Implementation of magneto-static problem (basic Implementation)
13 * \ingroup maxwell_element
14 *
15 * Look for theory and details here:
16 *
17 * \cite ivanyshyn2013computation
18 * <www.hpfem.jku.at/publications/szthesis.pdf>
19 * <https://pdfs.semanticscholar.org/0574/e5763d9b64bff16f908e9621f23af3b3dc86.pdf>
20 *
21 * Election file and all other problem related file are here \ref
22 * maxwell_element.
23 *
24 * \todo Extension for mix formulation
25 * \todo Use appropriate pre-conditioner for large problems
26 *
27 */
28struct MagneticElement {
29
30 MoFEM::Interface &mField;
31
32 /// \brief definition of volume element
33 struct VolumeFE : public MoFEM::VolumeElementForcesAndSourcesCore {
36 int getRule(int order) { return 2 * order + 1; };
37 };
38
39 // /// \brief definition of volume element
40 // struct VolumeFEReducedIntegration: public
41 // MoFEM::VolumeElementForcesAndSourcesCore {
42 // VolumeFEReducedIntegration(MoFEM::Interface &m_field):
43 // MoFEM::VolumeElementForcesAndSourcesCore(m_field) {}
44 // int getRule(int order) { return 2*order+1; };
45 // };
46
47 /** \brief define surface element
48 *
49 */
50 struct TriFE : public MoFEM::FaceElementForcesAndSourcesCore {
53 int getRule(int order) { return 2 * order; };
54 };
55
56 MagneticElement(MoFEM::Interface &m_field) : mField(m_field) {}
57 virtual ~MagneticElement() = default;
58
59 /**
60 * \brief data structure storing material constants, model parameters,
61 * matrices, etc.
62 *
63 */
64 struct BlockData {
65
66 // field
67 const string fieldName;
68 const string feName;
69 const string feNaturalBCName;
70
71 // material parameters
72 double mU; ///< magnetic constant N / A2
73 double ePsilon; ///< regularization paramater
74
75 // Natural boundary conditions
76 Range naturalBc;
77
78 // Essential boundary conditions
79 Range essentialBc;
80
81 int oRder; ///< approximation order
82
83 // Petsc data
84 DM dM;
85 Mat A;
86 Vec D, F;
87
88 BlockData()
89 : fieldName("MAGNETIC_POTENTIAL"), feName("MAGNETIC"),
90 feNaturalBCName("MAGENTIC_NATURAL_BC"), mU(1), ePsilon(0.1) {}
91 ~BlockData() {}
92 };
93
94 BlockData blockData;
95
96 /**
97 * \brief get natural boundary conditions
98 * \ingroup maxwell_element
99 * @return error code
100 */
101 MoFEMErrorCode getNaturalBc() {
102 MoFEMFunctionBegin;
103 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, bit)) {
104 if (bit->getName().compare(0, 9, "NATURALBC") == 0) {
105 Range faces;
106 CHKERR mField.get_moab().get_entities_by_type(bit->meshset, MBTRI,
107 faces, true);
108 CHKERR mField.get_moab().get_adjacencies(
109 faces, 1, true, blockData.naturalBc, moab::Interface::UNION);
110 blockData.naturalBc.merge(faces);
111 }
112 }
113 MoFEMFunctionReturn(0);
114 }
115
116 /**
117 * \brief get essential boundary conditions (only homogenous case is
118 * considered) \ingroup maxwell_element
119 * @return error code
120 */
121 MoFEMErrorCode getEssentialBc() {
122 ParallelComm *pcomm =
123 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
124 MoFEMFunctionBegin;
125 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, bit)) {
126 if (bit->getName().compare(0, 10, "ESSENTIALBC") == 0) {
127 Range faces;
128 CHKERR mField.get_moab().get_entities_by_type(bit->meshset, MBTRI,
129 faces, true);
130 CHKERR mField.get_moab().get_adjacencies(
131 faces, 1, true, blockData.essentialBc, moab::Interface::UNION);
132 blockData.essentialBc.merge(faces);
133 }
134 }
135 if (blockData.essentialBc.empty()) {
136 Range tets;
137 CHKERR mField.get_moab().get_entities_by_type(0, MBTET, tets);
138 Skinner skin(&mField.get_moab());
139 Range skin_faces; // skin faces from 3d ents
140 CHKERR skin.find_skin(0, tets, false, skin_faces);
141 skin_faces = subtract(skin_faces, blockData.naturalBc);
142 Range proc_skin;
143 CHKERR pcomm->filter_pstatus(skin_faces,
144 PSTATUS_SHARED | PSTATUS_MULTISHARED,
145 PSTATUS_NOT, -1, &proc_skin);
146 CHKERR mField.get_moab().get_adjacencies(
147 proc_skin, 1, true, blockData.essentialBc, moab::Interface::UNION);
148 blockData.essentialBc.merge(proc_skin);
149 }
150 MoFEMFunctionReturn(0);
151 }
152
153 /**
154 * \brief build problem data structures
155 * \ingroup maxwell_element
156 * @return error code
157 */
158 MoFEMErrorCode createFields() {
159 MoFEMFunctionBegin;
160
161 // Set entities bit level. each entity has bit level depending for example
162 // on refinement level. In this case we do not refine mesh or not do
163 // topological changes, simply set refinement level to zero on all entities.
164
165 CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(
166 0, 3, BitRefLevel().set(0));
167
168 // add fields
169 CHKERR mField.add_field(blockData.fieldName, HCURL, DEMKOWICZ_JACOBI_BASE,
170 1);
171 CHKERR mField.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
172 3);
173 // meshset consisting all entities in mesh
174 EntityHandle root_set = mField.get_moab().get_root_set();
175 // add entities to field
176 CHKERR mField.add_ents_to_field_by_type(root_set, MBTET,
177 blockData.fieldName);
178
179 // // The higher-order gradients can be gauged by locally skipping the
180 // // corresponding degrees of freedom and basis functions in the
181 // higher-order
182 // // edge-based, face-based and cell-based finite element subspaces.
183 //
184 // Range tris,edges;
185 // CHKERR mField.get_moab().get_entities_by_type(root_set,MBTRI,tris,true);
186 // mField.get_moab().get_entities_by_type(root_set,MBEDGE,edges,true);
187 //
188 // // Set order in volume
189 // Range bc_ents = unite(blockData.naturalBc,blockData.essentialBc);
190 // Range vol_ents = subtract(unite(tris,edges),bc_ents);
191 // CHKERR
192 // mField.set_field_order(vol_ents,blockData.fieldName,blockData.oRder); int
193 // gauged_order = 1; CHKERR
194 // mField.set_field_order(bc_ents,blockData.fieldName,gauged_order);
195
196 // Set order on tets
197 CHKERR mField.set_field_order(root_set, MBTET, blockData.fieldName,
198 blockData.oRder);
199 CHKERR mField.set_field_order(root_set, MBTRI, blockData.fieldName,
200 blockData.oRder);
201 CHKERR mField.set_field_order(root_set, MBEDGE, blockData.fieldName,
202 blockData.oRder);
203
204 // Set geometry approximation ordered
205 CHKERR mField.add_ents_to_field_by_type(root_set, MBTET,
206 "MESH_NODE_POSITIONS");
207 CHKERR mField.set_field_order(root_set, MBTET, "MESH_NODE_POSITIONS", 2);
208 CHKERR mField.set_field_order(root_set, MBTRI, "MESH_NODE_POSITIONS", 2);
209 CHKERR mField.set_field_order(root_set, MBEDGE, "MESH_NODE_POSITIONS", 2);
210 CHKERR mField.set_field_order(root_set, MBVERTEX, "MESH_NODE_POSITIONS", 1);
211
212 // build field
213 CHKERR mField.build_fields();
214
215 // get HO geometry for 10 node tets
216 // This method takes coordinates form edges mid nodes in 10 node tet and
217 // project values on 2nd order hierarchical basis used to approx. geometry.
218 Projection10NodeCoordsOnField ent_method_material(mField,
219 "MESH_NODE_POSITIONS");
220 CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
221
222 MoFEMFunctionReturn(0);
223 }
224
225 /**
226 * \brief create finite elements
227 * \ingroup maxwell_element
228 *
229 * Create volume and surface element. Surface element is used to integrate
230 * natural boundary conditions.
231 *
232 * @return error code
233 */
265
266 /**
267 * \brief create problem
268 * \ingroup maxwell_element
269 *
270 * Problem is collection of finite elements. With the information on which
271 * fields finite elements operates the matrix and left and right hand side
272 * vector could be created.
273 *
274 * Here we use Distributed mesh manager from PETSc as a inteface.
275 *
276 * @return error code
277 */
278 MoFEMErrorCode createProblem() {
279 MoFEMFunctionBegin;
280 // set up DM
281 CHKERR DMRegister_MoFEM("DMMOFEM");
282 CHKERR DMCreate(PETSC_COMM_WORLD, &blockData.dM);
283 CHKERR DMSetType(blockData.dM, "DMMOFEM");
284 CHKERR DMMoFEMCreateMoFEM(blockData.dM, &mField, "MAGNETIC_PROBLEM",
285 BitRefLevel().set(0));
286 CHKERR DMSetFromOptions(blockData.dM);
287 CHKERR DMMoFEMSetIsPartitioned(blockData.dM, PETSC_TRUE);
288 // add elements to blockData.dM
289 CHKERR DMMoFEMAddElement(blockData.dM, blockData.feName);
290 CHKERR DMMoFEMAddElement(blockData.dM, blockData.feNaturalBCName);
291 CHKERR DMSetUp(blockData.dM);
292
293 // remove essential DOFs
294 const MoFEM::Problem *problem_ptr;
295 CHKERR DMMoFEMGetProblemPtr(blockData.dM, &problem_ptr);
296 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
297 problem_ptr->getName(), blockData.fieldName, blockData.essentialBc);
298
299 MoFEMFunctionReturn(0);
300 }
301
302 /**
303 * \brief destroy Distributed mesh manager
304 * \ingroup maxwell_element
305 * @return [description]
306 */
307 MoFEMErrorCode destroyProblem() {
308 MoFEMFunctionBegin;
309 CHKERR DMDestroy(&blockData.dM);
310 MoFEMFunctionReturn(0);
311 }
312
313 /** \brief solve problem
314 * \ingroup maxwell_element
315 *
316 * Create matrices; integrate over elements; solve linear system of equations
317 *
318 */
319 MoFEMErrorCode solveProblem(const bool regression_test = false) {
320 MoFEMFunctionBegin;
321
322 VolumeFE vol_fe(mField);
323 auto material_grad_mat = boost::make_shared<MatrixDouble>();
324 auto material_det_vec = boost::make_shared<VectorDouble>();
325 auto material_inv_grad_mat = boost::make_shared<MatrixDouble>();
326 CHKERR AddHOOps<3, 3, 3>::add(vol_fe.getOpPtrVector(), {HCURL, L2},
327 "MESH_NODE_POSITIONS");
328 vol_fe.getOpPtrVector().push_back(new OpCurlCurl(blockData));
329 vol_fe.getOpPtrVector().push_back(new OpStab(blockData));
330 TriFE tri_fe(mField);
331 tri_fe.meshPositionsFieldName = "none";
332
333 tri_fe.getOpPtrVector().push_back(
334 new OpGetHONormalsOnFace("MESH_NODE_POSITIONS"));
335 tri_fe.getOpPtrVector().push_back(
336 new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
337 tri_fe.getOpPtrVector().push_back(
338 new OpHOSetCovariantPiolaTransformOnFace3D(HCURL));
339 tri_fe.getOpPtrVector().push_back(new OpNaturalBC(blockData));
340
341 // create matrices and vectors
342 CHKERR DMCreateGlobalVector(blockData.dM, &blockData.D);
343 // CHKERR DMCreateMatrix(blockData.dM, &blockData.A);
344 CHKERR DMCreateGlobalVector(blockData.dM, &blockData.F);
345 CHKERR VecDuplicate(blockData.F, &blockData.D);
346
347 const MoFEM::Problem *problem_ptr;
348 CHKERR DMMoFEMGetProblemPtr(blockData.dM, &problem_ptr);
349 CHKERR mField.getInterface<MatrixManager>()
350 ->createMPIAIJ<PetscGlobalIdx_mi_tag>(problem_ptr->getName(),
351 &blockData.A);
352
353 CHKERR MatZeroEntries(blockData.A);
354 CHKERR VecZeroEntries(blockData.F);
355 CHKERR VecGhostUpdateBegin(blockData.F, ADD_VALUES, SCATTER_FORWARD);
356 CHKERR VecGhostUpdateEnd(blockData.F, ADD_VALUES, SCATTER_FORWARD);
357
358 CHKERR DMoFEMLoopFiniteElements(blockData.dM, blockData.feName.c_str(),
359 &vol_fe);
360 CHKERR DMoFEMLoopFiniteElements(blockData.dM,
361 blockData.feNaturalBCName.c_str(), &tri_fe);
362
363 CHKERR MatAssemblyBegin(blockData.A, MAT_FINAL_ASSEMBLY);
364 CHKERR MatAssemblyEnd(blockData.A, MAT_FINAL_ASSEMBLY);
365 CHKERR VecGhostUpdateBegin(blockData.F, ADD_VALUES, SCATTER_REVERSE);
366 CHKERR VecGhostUpdateEnd(blockData.F, ADD_VALUES, SCATTER_REVERSE);
367 CHKERR VecAssemblyBegin(blockData.F);
368 CHKERR VecAssemblyEnd(blockData.F);
369
370 KSP solver;
371 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
372 CHKERR KSPSetOperators(solver, blockData.A, blockData.A);
373 CHKERR KSPSetFromOptions(solver);
374 CHKERR KSPSetUp(solver);
375 CHKERR KSPSolve(solver, blockData.F, blockData.D);
376 CHKERR KSPDestroy(&solver);
377
378 CHKERR VecGhostUpdateBegin(blockData.D, INSERT_VALUES, SCATTER_FORWARD);
379 CHKERR VecGhostUpdateEnd(blockData.D, INSERT_VALUES, SCATTER_FORWARD);
380 CHKERR DMoFEMMeshToLocalVector(blockData.dM, blockData.D, INSERT_VALUES,
381 SCATTER_REVERSE);
382
383 if (regression_test) {
384 // This test is for order = 1 only
385 double nrm2;
386 CHKERR VecNorm(blockData.D, NORM_2, &nrm2);
387 const double expected_value = 4.6772e+01;
388 const double tol = 1e-4;
389 if ((nrm2 - expected_value) / expected_value > tol) {
390 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
391 "Regression test failed %6.4e != %6.4e", nrm2, expected_value);
392 }
393 }
394
395 CHKERR VecDestroy(&blockData.D);
396 CHKERR VecDestroy(&blockData.F);
397 CHKERR MatDestroy(&blockData.A);
398
399 MoFEMFunctionReturn(0);
400 }
401
402 /**
403 * \brief post-process results, i.e. save solution on the mesh
404 * \ingroup maxwell_element
405 * @return [description]
406 */
407 MoFEMErrorCode postProcessResults() {
408
409 MoFEMFunctionBegin;
410 PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore> post_proc(
411 mField);
412
413 CHKERR AddHOOps<3, 3, 3>::add(post_proc.getOpPtrVector(), {HCURL, L2},
414 "MESH_NODE_POSITIONS");
415
416 auto pos_ptr = boost::make_shared<MatrixDouble>();
417 auto field_val_ptr = boost::make_shared<MatrixDouble>();
418
419 post_proc.getOpPtrVector().push_back(
420 new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS", pos_ptr));
421 post_proc.getOpPtrVector().push_back(
422 new OpCalculateHVecVectorField<3>(blockData.fieldName, field_val_ptr));
423
424 using OpPPMap = OpPostProcMapInMoab<3, 3>;
425
426 post_proc.getOpPtrVector().push_back(
427
428 new OpPPMap(
429
430 post_proc.getPostProcMesh(), post_proc.getMapGaussPts(),
431
432 {},
433
434 {{"MESH_NODE_POSITIONS", pos_ptr},
435 {blockData.fieldName, field_val_ptr}},
436
437 {},
438
439 {}
440
441 )
442
443 );
444
445 post_proc.getOpPtrVector().push_back(new OpPostProcessCurl(
446 blockData, post_proc.getPostProcMesh(), post_proc.getMapGaussPts()));
447 CHKERR DMoFEMLoopFiniteElements(blockData.dM, blockData.feName.c_str(),
448 &post_proc);
449 CHKERR post_proc.writeFile("out_values.h5m");
450 MoFEMFunctionReturn(0);
451 }
452
453 /** \brief calculate and assemble CurlCurl matrix
454 * \ingroup maxwell_element
455
456 \f[
457 \mathbf{A} = \int_\Omega \mu^{-1} \left( \nabla \times \mathbf{u} \cdot
458 \nabla \times \mathbf{v} \right) \textrm{d}\Omega \f] where \f[ \mathbf{u} =
459 \nabla \times \mathbf{B} \f] where \f$\mathbf{B}\f$ is magnetic flux and
460 \f$\mu\f$ is magnetic permeability.
461
462 For more details pleas look to \cite ivanyshyn2013computation
463
464 */
465 struct OpCurlCurl
466 : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
467
468 BlockData &blockData;
469 OpCurlCurl(BlockData &data)
470 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
471 data.fieldName, UserDataOperator::OPROWCOL),
472 blockData(data) {
473 sYmm = true;
474 }
475
476 MatrixDouble entityLocalMatrix;
477
478 /**
479 * \brief integrate matrix
480 * @param row_side local number of entity on element for row of the matrix
481 * @param col_side local number of entity on element for col of the matrix
482 * @param row_type type of row entity (EDGE/TRIANGLE/TETRAHEDRON)
483 * @param col_type type of col entity (EDGE/TRIANGLE/TETRAHEDRON)
484 * @param row_data structure of data, like base functions and associated
485 * methods to access those data on rows
486 * @param col_data structure of data, like base functions and associated
487 * methods to access those data on rows
488 * @return error code
489 */
490 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
491 EntityType col_type,
492 EntitiesFieldData::EntData &row_data,
493 EntitiesFieldData::EntData &col_data) {
494 MoFEMFunctionBegin;
495
496 if (row_type == MBVERTEX)
497 MoFEMFunctionReturnHot(0);
498 if (col_type == MBVERTEX)
499 MoFEMFunctionReturnHot(0);
500
501 const int nb_row_dofs = row_data.getN().size2() / 3;
502 if (nb_row_dofs == 0)
503 MoFEMFunctionReturnHot(0);
504 const int nb_col_dofs = col_data.getN().size2() / 3;
505 if (nb_col_dofs == 0)
506 MoFEMFunctionReturnHot(0);
507 entityLocalMatrix.resize(nb_row_dofs, nb_col_dofs, false);
508 entityLocalMatrix.clear();
509
510 if (nb_row_dofs != static_cast<int>(row_data.getFieldData().size())) {
511 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
512 "Number of base functions and DOFs on entity is different on "
513 "rows %d!=%zu",
514 nb_row_dofs, row_data.getFieldData().size());
515 }
516 if (nb_col_dofs != static_cast<int>(col_data.getFieldData().size())) {
517 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
518 "Number of base functions and DOFs on entity is different on "
519 "cols %d!=%zu",
520 nb_col_dofs, col_data.getFieldData().size());
521 }
522
523 MatrixDouble row_curl_mat, col_curl_mat;
524 FTensor::Index<'i', 3> i;
525 FTensor::Index<'j', 3> j;
526 FTensor::Index<'k', 3> k;
527
528 const double c0 = 1. / blockData.mU;
529 const int nb_gauss_pts = row_data.getN().size1();
530 auto t_row_curl_base = row_data.getFTensor2DiffN<3, 3>();
531
532 for (int gg = 0; gg != nb_gauss_pts; gg++) {
533
534 // get integration weight scaled by volume
535 double w = getGaussPts()(3, gg) * getVolume();
536
537 FTensor::Tensor1<double, 3> t_row_curl;
538 for (int aa = 0; aa != nb_row_dofs; aa++) {
539 t_row_curl(i) = levi_civita(j, i, k) * t_row_curl_base(j, k);
540
541 FTensor::Tensor0<double *> t_local_mat(&entityLocalMatrix(aa, 0), 1);
542 FTensor::Tensor1<double, 3> t_col_curl;
543
544 auto t_col_curl_base = col_data.getFTensor2DiffN<3, 3>(gg, 0);
545 for (int bb = 0; bb != nb_col_dofs; bb++) {
546 t_col_curl(i) = levi_civita(j, i, k) * t_col_curl_base(j, k);
547 t_local_mat += c0 * w * t_row_curl(i) * t_col_curl(i);
548 ++t_local_mat;
549 ++t_col_curl_base;
550 }
551
552 ++t_row_curl_base;
553 }
554 }
555
556 CHKERR MatSetValues(blockData.A, nb_row_dofs, &row_data.getIndices()[0],
557 nb_col_dofs, &col_data.getIndices()[0],
558 &entityLocalMatrix(0, 0), ADD_VALUES);
559
560 if (row_side != col_side || row_type != col_type) {
561 entityLocalMatrix = trans(entityLocalMatrix);
562 CHKERR MatSetValues(blockData.A, nb_col_dofs, &col_data.getIndices()[0],
563 nb_row_dofs, &row_data.getIndices()[0],
564 &entityLocalMatrix(0, 0), ADD_VALUES);
565 }
566
567 MoFEMFunctionReturn(0);
568 }
569 };
570
571 /** \brief calculate and assemble stabilization matrix
572 * \ingroup maxwell_element
573
574 \f[
575 \mathbf{A} = \int_\Omega \epsilon \mathbf{u} \cdot \mathbf{v}
576 \textrm{d}\Omega \f] where \f$\epsilon\f$ is regularization parameter.
577
578 */
579 struct OpStab
580 : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
581
582 BlockData &blockData;
583 OpStab(BlockData &data)
584 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
585 data.fieldName, UserDataOperator::OPROWCOL),
586 blockData(data) {
587 sYmm = true;
588 }
589
590 MatrixDouble entityLocalMatrix;
591
592 /**
593 * \brief integrate matrix
594 * @param row_side local number of entity on element for row of the matrix
595 * @param col_side local number of entity on element for col of the matrix
596 * @param row_type type of row entity (EDGE/TRIANGLE/TETRAHEDRON)
597 * @param col_type type of col entity (EDGE/TRIANGLE/TETRAHEDRON)
598 * @param row_data structure of data, like base functions and associated
599 * methods to access those data on rows
600 * @param col_data structure of data, like base functions and associated
601 * methods to access those data on rows
602 * @return error code
603 */
604 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
605 EntityType col_type,
606 EntitiesFieldData::EntData &row_data,
607 EntitiesFieldData::EntData &col_data) {
608 MoFEMFunctionBegin;
609
610 if (row_type == MBVERTEX)
611 MoFEMFunctionReturnHot(0);
612 if (col_type == MBVERTEX)
613 MoFEMFunctionReturnHot(0);
614
615 const int nb_row_dofs = row_data.getN().size2() / 3;
616 if (nb_row_dofs == 0)
617 MoFEMFunctionReturnHot(0);
618 const int nb_col_dofs = col_data.getN().size2() / 3;
619 if (nb_col_dofs == 0)
620 MoFEMFunctionReturnHot(0);
621 entityLocalMatrix.resize(nb_row_dofs, nb_col_dofs, false);
622 entityLocalMatrix.clear();
623
624 if (nb_row_dofs != static_cast<int>(row_data.getFieldData().size())) {
625 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
626 "Number of base functions and DOFs on entity is different on "
627 "rows %d!=%zu",
628 nb_row_dofs, row_data.getFieldData().size());
629 }
630 if (nb_col_dofs != static_cast<int>(col_data.getFieldData().size())) {
631 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
632 "Number of base functions and DOFs on entity is different on "
633 "cols %d!=%zu",
634 nb_col_dofs, col_data.getFieldData().size());
635 }
636
637 MatrixDouble row_curl_mat, col_curl_mat;
638 FTensor::Index<'i', 3> i;
639
640 const double c0 = 1. / blockData.mU;
641 const double c1 = blockData.ePsilon * c0;
642 const int nb_gauss_pts = row_data.getN().size1();
643
644 for (int gg = 0; gg != nb_gauss_pts; gg++) {
645
646 // get integration weight scaled by volume
647 double w = getGaussPts()(3, gg) * getVolume();
648
649 FTensor::Tensor1<const double *, 3> t_row_base(
650 &row_data.getVectorN<3>(gg)(0, HVEC0),
651 &row_data.getVectorN<3>(gg)(0, HVEC1),
652 &row_data.getVectorN<3>(gg)(0, HVEC2), 3);
653
654 for (int aa = 0; aa != nb_row_dofs; aa++) {
655 FTensor::Tensor0<double *> t_local_mat(&entityLocalMatrix(aa, 0), 1);
656 FTensor::Tensor1<const double *, 3> t_col_base(
657 &col_data.getVectorN<3>(gg)(0, HVEC0),
658 &col_data.getVectorN<3>(gg)(0, HVEC1),
659 &col_data.getVectorN<3>(gg)(0, HVEC2), 3);
660 for (int bb = 0; bb != nb_col_dofs; bb++) {
661 t_local_mat += c1 * w * t_row_base(i) * t_col_base(i);
662 ++t_col_base;
663 ++t_local_mat;
664 }
665 ++t_row_base;
666 }
667 }
668
669 // cerr << entityLocalMatrix << endl;
670 // cerr << endl;
671
672 CHKERR MatSetValues(blockData.A, nb_row_dofs, &row_data.getIndices()[0],
673 nb_col_dofs, &col_data.getIndices()[0],
674 &entityLocalMatrix(0, 0), ADD_VALUES);
675
676 if (row_side != col_side || row_type != col_type) {
677 entityLocalMatrix = trans(entityLocalMatrix);
678 CHKERR MatSetValues(blockData.A, nb_col_dofs, &col_data.getIndices()[0],
679 nb_row_dofs, &row_data.getIndices()[0],
680 &entityLocalMatrix(0, 0), ADD_VALUES);
681 }
682
683 MoFEMFunctionReturn(0);
684 }
685 };
686
687 /** \brief calculate essential boundary conditions
688 * \ingroup maxwell_element
689
690 \f[
691 \mathbf{F} = \int_{\partial\Omega} \ \mathbf{u} \cdot \mathbf{j}_i
692 \textrm{d}{\partial\Omega} \f] where \f$\mathbf{j}_i\f$ is current density
693 function.
694
695 Here simple current on coil is hard coded. In future more general
696 implementation is needed.
697
698 */
699 struct OpNaturalBC
700 : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
701
702 BlockData &blockData;
703
704 OpNaturalBC(BlockData &data)
705 : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
706 data.fieldName, UserDataOperator::OPROW),
707 blockData(data) {}
708
709 VectorDouble naturalBC;
710
711 /**
712 * \brief integrate matrix
713 * \ingroup maxwell_element
714 * @param row_side local number of entity on element for row of the matrix
715 * @param row_type type of row entity (EDGE/TRIANGLE/TETRAHEDRON)
716 * @param row_data structure of data, like base functions and associated
717 * methods to access those data on rows
718 * @return error code
719 */
720 MoFEMErrorCode doWork(int row_side, EntityType row_type,
721 EntitiesFieldData::EntData &row_data) {
722 MoFEMFunctionBegin;
723
724 if (row_type == MBVERTEX)
725 MoFEMFunctionReturnHot(0);
726
727 const int nb_row_dofs = row_data.getN().size2() / 3;
728 if (nb_row_dofs == 0)
729 MoFEMFunctionReturnHot(0);
730 naturalBC.resize(nb_row_dofs, false);
731 naturalBC.clear();
732
733 FTensor::Index<'i', 3> i;
734
735 const int nb_gauss_pts = row_data.getN().size1();
736 auto t_row_base = row_data.getFTensor1N<3>();
737
738 for (int gg = 0; gg != nb_gauss_pts; gg++) {
739
740 // get integration weight scaled by volume
741 double area;
742 area = norm_2(getNormalsAtGaussPts(gg)) * 0.5;
743 double w = getGaussPts()(2, gg) * area;
744
745 // Current is on surface where natural bc are applied. It is set that
746 // current is in XY plane, circular, around the coil.
747 const double x = getCoordsAtGaussPts()(gg, 0);
748 const double y = getCoordsAtGaussPts()(gg, 1);
749 const double r = sqrt(x * x + y * y);
750 FTensor::Tensor1<double, 3> t_j;
751 t_j(0) = -y / r;
752 t_j(1) = +x / r;
753 t_j(2) = 0;
754
755 // double a = t_j(i) * t_tangent1(i);
756 // double b = t_j(i) * t_tangent2(i);
757 // t_j(i) = a * t_tangent1(i) + b * t_tangent2(i);
758
759 // ++t_tangent1;
760 // ++t_tangent2;
761
762 FTensor::Tensor0<double *> t_f(&naturalBC[0]);
763 for (int aa = 0; aa != nb_row_dofs; aa++) {
764 t_f += w * t_row_base(i) * t_j(i);
765 ++t_row_base;
766 ++t_f;
767 }
768 }
769
770 CHKERR VecSetValues(blockData.F, row_data.getIndices().size(),
771 &row_data.getIndices()[0], &naturalBC[0], ADD_VALUES);
772
773 MoFEMFunctionReturn(0);
774 }
775 };
776
777 /** \brief calculate and assemble CurlCurl matrix
778 * \ingroup maxwell_element
779 */
780 struct OpPostProcessCurl
781 : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
782
783 BlockData &blockData;
784 moab::Interface &postProcMesh;
785 std::vector<EntityHandle> &mapGaussPts;
786
787 OpPostProcessCurl(BlockData &data, moab::Interface &post_proc_mesh,
788 std::vector<EntityHandle> &map_gauss_pts)
789 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
790 data.fieldName, UserDataOperator::OPROW),
791 blockData(data), postProcMesh(post_proc_mesh),
792 mapGaussPts(map_gauss_pts) {}
793
794 MoFEMErrorCode doWork(int row_side, EntityType row_type,
795 EntitiesFieldData::EntData &row_data) {
796 MoFEMFunctionBegin;
797
798 if (row_type == MBVERTEX)
799 MoFEMFunctionReturnHot(0);
800
801 Tag th;
802 double def_val[] = {0, 0, 0};
803 CHKERR postProcMesh.tag_get_handle("MAGNETIC_INDUCTION_FIELD", 3,
804 MB_TYPE_DOUBLE, th,
805 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
806 const int nb_row_dofs = row_data.getN().size2() / 3;
807 if (nb_row_dofs == 0)
808 MoFEMFunctionReturnHot(0);
809 const void *tags_ptr[mapGaussPts.size()];
810
811 if (nb_row_dofs != row_data.getFieldData().size())
812 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
813 "Wrong number of base functions and DOFs %d != %zu",
814 nb_row_dofs, row_data.getFieldData().size());
815
816 FTensor::Index<'i', 3> i;
817 FTensor::Index<'j', 3> j;
818 FTensor::Index<'k', 3> k;
819 const int nb_gauss_pts = row_data.getN().size1();
820 if (nb_gauss_pts != static_cast<int>(mapGaussPts.size())) {
821 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
822 "Inconsistency number of dofs %d!=%zu", nb_gauss_pts,
823 mapGaussPts.size());
824 }
825 CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0],
826 mapGaussPts.size(), tags_ptr);
827 auto t_curl_base = row_data.getFTensor2DiffN<3, 3>();
828 for (int gg = 0; gg != nb_gauss_pts; gg++) {
829 // get pointer to tag values on entity (i.e. vertex on refined
830 // post-processing mesh)
831 double *ptr = &((double *)tags_ptr[gg])[0];
832 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_curl(ptr, &ptr[1],
833 &ptr[2]);
834 // calculate curl value
835 auto t_dof = row_data.getFTensor0FieldData();
836 for (int aa = 0; aa != nb_row_dofs; aa++) {
837 t_curl(i) += t_dof * (levi_civita(j, i, k) * t_curl_base(j, k));
838 ++t_curl_base;
839 ++t_dof;
840 }
841 }
842 MoFEMFunctionReturn(0);
843 }
844 };
845};
846
847#endif //__MAGNETICELEMENT_HPP__
848
849/******************************************************************************
850 * \defgroup maxwell_element Magnetic/Maxwell element
851 * \ingroup user_modules
852 ******************************************************************************/
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition HookeElement.hpp:75
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
H1
@ H1
continuous field
Definition definitions.h:85
HCURL
@ HCURL
field with continuous tangents
Definition definitions.h:86
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition definitions.h:228
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
BLOCKSET
@ BLOCKSET
Definition definitions.h:161
HVEC0
@ HVEC0
Definition definitions.h:199
HVEC1
@ HVEC1
Definition definitions.h:199
HVEC2
@ HVEC2
Definition definitions.h:199
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
order
constexpr int order
Definition dg_projection.cpp:38
MoFEM::DMoFEMLoopFiniteElements
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
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::CoreInterface::modify_finite_element_add_field_col
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string name, const bool recursive=true)=0
add entities to finite element
MoFEM::CoreInterface::modify_finite_element_add_field_row
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::CoreInterface::set_field_order
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
MoFEM::CoreInterface::add_ents_to_field_by_type
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
MoFEM::CoreInterface::loop_dofs
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.
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#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.
Definition MeshsetsManager.hpp:94
bit
auto bit
set bit
Definition hanging_node_approx.cpp:75
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:21
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:22
tol
double tol
Definition mesh_smoothing.cpp:23
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition photon_diffusion.cpp:40
MagneticElement::BlockData
data structure storing material constants, model parameters, matrices, etc.
Definition MagneticElement.hpp:64
MagneticElement::BlockData::oRder
int oRder
approximation order
Definition MagneticElement.hpp:81
MagneticElement::BlockData::mU
double mU
magnetic constant N / A2
Definition MagneticElement.hpp:72
MagneticElement::BlockData::ePsilon
double ePsilon
regularization paramater
Definition MagneticElement.hpp:73
MagneticElement::OpCurlCurl
calculate and assemble CurlCurl matrix
Definition MagneticElement.hpp:466
MagneticElement::OpCurlCurl::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
integrate matrix
Definition MagneticElement.hpp:490
MagneticElement::OpNaturalBC
calculate essential boundary conditions
Definition MagneticElement.hpp:700
MagneticElement::OpNaturalBC::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
integrate matrix
Definition MagneticElement.hpp:720
MagneticElement::OpPostProcessCurl
calculate and assemble CurlCurl matrix
Definition MagneticElement.hpp:781
MagneticElement::OpPostProcessCurl::mapGaussPts
std::vector< EntityHandle > & mapGaussPts
Definition MagneticElement.hpp:785
MagneticElement::OpPostProcessCurl::OpPostProcessCurl
OpPostProcessCurl(BlockData &data, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts)
Definition MagneticElement.hpp:787
MagneticElement::OpPostProcessCurl::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Definition MagneticElement.hpp:794
MagneticElement::OpStab
calculate and assemble stabilization matrix
Definition MagneticElement.hpp:580
MagneticElement::OpStab::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
integrate matrix
Definition MagneticElement.hpp:604
MagneticElement::TriFE
define surface element
Definition MagneticElement.hpp:50
MagneticElement::TriFE::TriFE
TriFE(MoFEM::Interface &m_field)
Definition MagneticElement.hpp:51
MagneticElement::VolumeFE
definition of volume element
Definition MagneticElement.hpp:33
MagneticElement::VolumeFE::VolumeFE
VolumeFE(MoFEM::Interface &m_field)
Definition MagneticElement.hpp:34
MagneticElement
Implementation of magneto-static problem (basic Implementation)
Definition MagneticElement.hpp:28
MagneticElement::createFields
MoFEMErrorCode createFields()
build problem data structures
Definition MagneticElement.hpp:158
MagneticElement::createElements
MoFEMErrorCode createElements()
create finite elements
Definition MagneticElement.hpp:234
MagneticElement::destroyProblem
MoFEMErrorCode destroyProblem()
destroy Distributed mesh manager
Definition MagneticElement.hpp:307
MagneticElement::postProcessResults
MoFEMErrorCode postProcessResults()
post-process results, i.e. save solution on the mesh
Definition MagneticElement.hpp:407
MagneticElement::createProblem
MoFEMErrorCode createProblem()
create problem
Definition MagneticElement.hpp:278
MagneticElement::mField
MoFEM::Interface & mField
Definition MagneticElement.hpp:30
MagneticElement::getEssentialBc
MoFEMErrorCode getEssentialBc()
get essential boundary conditions (only homogenous case is considered)
Definition MagneticElement.hpp:121
MagneticElement::getNaturalBc
MoFEMErrorCode getNaturalBc()
get natural boundary conditions
Definition MagneticElement.hpp:101
MagneticElement::MagneticElement
MagneticElement(MoFEM::Interface &m_field)
Definition MagneticElement.hpp:56
MagneticElement::~MagneticElement
virtual ~MagneticElement()=default
MagneticElement::solveProblem
MoFEMErrorCode solveProblem(const bool regression_test=false)
solve problem
Definition MagneticElement.hpp:319
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEM::CoreInterface::add_field
virtual MoFEMErrorCode add_field(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_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition ForcesAndSourcesCore.hpp:121
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition PostProcBrokenMeshInMoabBase.hpp:745
MoFEM::Problem
keeps basic data about problem
Definition ProblemsMultiIndices.hpp:54
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93
Generated by Doxygen 1.9.8 and hosted at University of Glasgow