v0.14.0
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 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  */
28 struct MagneticElement {
29 
30  MoFEM::Interface &mField;
31 
32  /// \brief definition of volume element
33  struct VolumeFE : public MoFEM::VolumeElementForcesAndSourcesCore {
34  VolumeFE(MoFEM::Interface &m_field)
35  : MoFEM::VolumeElementForcesAndSourcesCore(m_field) {}
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 {
51  TriFE(MoFEM::Interface &m_field)
52  : MoFEM::FaceElementForcesAndSourcesCore(m_field) {}
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  */
234  MoFEMErrorCode createElements() {
235  MoFEMFunctionBegin;
236  // //Elements
237  CHKERR mField.add_finite_element(blockData.feName);
238  CHKERR mField.modify_finite_element_add_field_row(blockData.feName,
239  blockData.fieldName);
240  CHKERR mField.modify_finite_element_add_field_col(blockData.feName,
241  blockData.fieldName);
242  CHKERR mField.modify_finite_element_add_field_data(blockData.feName,
243  blockData.fieldName);
244  CHKERR mField.modify_finite_element_add_field_data(blockData.feName,
245  "MESH_NODE_POSITIONS");
246  CHKERR mField.add_finite_element(blockData.feNaturalBCName);
247  CHKERR mField.modify_finite_element_add_field_row(blockData.feNaturalBCName,
248  blockData.fieldName);
249  CHKERR mField.modify_finite_element_add_field_col(blockData.feNaturalBCName,
250  blockData.fieldName);
251  CHKERR mField.modify_finite_element_add_field_data(
252  blockData.feNaturalBCName, blockData.fieldName);
253  CHKERR mField.modify_finite_element_add_field_data(
254  blockData.feNaturalBCName, "MESH_NODE_POSITIONS");
255  CHKERR mField.add_ents_to_finite_element_by_type(0, MBTET,
256  blockData.feName);
257  CHKERR mField.add_ents_to_finite_element_by_type(blockData.naturalBc, MBTRI,
258  blockData.feNaturalBCName);
259  // build finite elemnts
260  CHKERR mField.build_finite_elements();
261  // build adjacencies
262  CHKERR mField.build_adjacencies(BitRefLevel().set(0));
263  MoFEMFunctionReturn(0);
264  }
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 addHOOpsVol("MESH_NODE_POSITIONS", vol_fe, false, true, false, true);
327  vol_fe.getOpPtrVector().push_back(new OpCurlCurl(blockData));
328  vol_fe.getOpPtrVector().push_back(new OpStab(blockData));
329  TriFE tri_fe(mField);
330  tri_fe.meshPositionsFieldName = "none";
331 
332  tri_fe.getOpPtrVector().push_back(
333  new OpGetHONormalsOnFace("MESH_NODE_POSITIONS"));
334  tri_fe.getOpPtrVector().push_back(
335  new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
336  tri_fe.getOpPtrVector().push_back(
337  new OpHOSetCovariantPiolaTransformOnFace3D(HCURL));
338  tri_fe.getOpPtrVector().push_back(new OpNaturalBC(blockData));
339 
340  // create matrices and vectors
341  CHKERR DMCreateGlobalVector(blockData.dM, &blockData.D);
342  // CHKERR DMCreateMatrix(blockData.dM, &blockData.A);
343  CHKERR DMCreateGlobalVector(blockData.dM, &blockData.F);
344  CHKERR VecDuplicate(blockData.F, &blockData.D);
345 
346  const MoFEM::Problem *problem_ptr;
347  CHKERR DMMoFEMGetProblemPtr(blockData.dM, &problem_ptr);
348  CHKERR mField.getInterface<MatrixManager>()
349  ->createMPIAIJ<PetscGlobalIdx_mi_tag>(problem_ptr->getName(),
350  &blockData.A);
351 
352  CHKERR MatZeroEntries(blockData.A);
353  CHKERR VecZeroEntries(blockData.F);
354  CHKERR VecGhostUpdateBegin(blockData.F, ADD_VALUES, SCATTER_FORWARD);
355  CHKERR VecGhostUpdateEnd(blockData.F, ADD_VALUES, SCATTER_FORWARD);
356 
357  CHKERR DMoFEMLoopFiniteElements(blockData.dM, blockData.feName.c_str(),
358  &vol_fe);
359  CHKERR DMoFEMLoopFiniteElements(blockData.dM,
360  blockData.feNaturalBCName.c_str(), &tri_fe);
361 
362  CHKERR MatAssemblyBegin(blockData.A, MAT_FINAL_ASSEMBLY);
363  CHKERR MatAssemblyEnd(blockData.A, MAT_FINAL_ASSEMBLY);
364  CHKERR VecGhostUpdateBegin(blockData.F, ADD_VALUES, SCATTER_REVERSE);
365  CHKERR VecGhostUpdateEnd(blockData.F, ADD_VALUES, SCATTER_REVERSE);
366  CHKERR VecAssemblyBegin(blockData.F);
367  CHKERR VecAssemblyEnd(blockData.F);
368 
369  KSP solver;
370  CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
371  CHKERR KSPSetOperators(solver, blockData.A, blockData.A);
372  CHKERR KSPSetFromOptions(solver);
373  CHKERR KSPSetUp(solver);
374  CHKERR KSPSolve(solver, blockData.F, blockData.D);
375  CHKERR KSPDestroy(&solver);
376 
377  CHKERR VecGhostUpdateBegin(blockData.D, INSERT_VALUES, SCATTER_FORWARD);
378  CHKERR VecGhostUpdateEnd(blockData.D, INSERT_VALUES, SCATTER_FORWARD);
379  CHKERR DMoFEMMeshToLocalVector(blockData.dM, blockData.D, INSERT_VALUES,
380  SCATTER_REVERSE);
381 
382  if (regression_test) {
383  // This test is for order = 1 only
384  double nrm2;
385  CHKERR VecNorm(blockData.D, NORM_2, &nrm2);
386  const double expected_value = 4.6772e+01;
387  const double tol = 1e-4;
388  if ((nrm2 - expected_value) / expected_value > tol) {
389  SETERRQ2(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
390  "Regression test failed %6.4e != %6.4e", nrm2, expected_value);
391  }
392  }
393 
394  CHKERR VecDestroy(&blockData.D);
395  CHKERR VecDestroy(&blockData.F);
396  CHKERR MatDestroy(&blockData.A);
397 
398  MoFEMFunctionReturn(0);
399  }
400 
401  /**
402  * \brief post-process results, i.e. save solution on the mesh
403  * \ingroup maxwell_element
404  * @return [description]
405  */
406  MoFEMErrorCode postProcessResults() {
407 
408  MoFEMFunctionBegin;
409  PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore> post_proc(
410  mField);
411 
412  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", post_proc, false, true, false,
413  true);
414 
415  auto pos_ptr = boost::make_shared<MatrixDouble>();
416  auto field_val_ptr = boost::make_shared<MatrixDouble>();
417 
418  post_proc.getOpPtrVector().push_back(
419  new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS", pos_ptr));
420  post_proc.getOpPtrVector().push_back(
421  new OpCalculateHVecVectorField<3>(blockData.fieldName, field_val_ptr));
422 
423  using OpPPMap = OpPostProcMapInMoab<3, 3>;
424 
425  post_proc.getOpPtrVector().push_back(
426 
427  new OpPPMap(
428 
429  post_proc.getPostProcMesh(), post_proc.getMapGaussPts(),
430 
431  {},
432 
433  {{"MESH_NODE_POSITIONS", pos_ptr},
434  {blockData.fieldName, field_val_ptr}},
435 
436  {},
437 
438  {}
439 
440  )
441 
442  );
443 
444  post_proc.getOpPtrVector().push_back(new OpPostProcessCurl(
445  blockData, post_proc.getPostProcMesh(), post_proc.getMapGaussPts()));
446  CHKERR DMoFEMLoopFiniteElements(blockData.dM, blockData.feName.c_str(),
447  &post_proc);
448  CHKERR post_proc.writeFile("out_values.h5m");
449  MoFEMFunctionReturn(0);
450  }
451 
452  /** \brief calculate and assemble CurlCurl matrix
453  * \ingroup maxwell_element
454 
455  \f[
456  \mathbf{A} = \int_\Omega \mu^{-1} \left( \nabla \times \mathbf{u} \cdot
457  \nabla \times \mathbf{v} \right) \textrm{d}\Omega \f] where \f[ \mathbf{u} =
458  \nabla \times \mathbf{B} \f] where \f$\mathbf{B}\f$ is magnetic flux and
459  \f$\mu\f$ is magnetic permeability.
460 
461  For more details pleas look to \cite ivanyshyn2013computation
462 
463  */
464  struct OpCurlCurl
465  : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
466 
467  BlockData &blockData;
468  OpCurlCurl(BlockData &data)
469  : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
470  data.fieldName, UserDataOperator::OPROWCOL),
471  blockData(data) {
472  sYmm = true;
473  }
474 
475  MatrixDouble entityLocalMatrix;
476 
477  /**
478  * \brief integrate matrix
479  * @param row_side local number of entity on element for row of the matrix
480  * @param col_side local number of entity on element for col of the matrix
481  * @param row_type type of row entity (EDGE/TRIANGLE/TETRAHEDRON)
482  * @param col_type type of col entity (EDGE/TRIANGLE/TETRAHEDRON)
483  * @param row_data structure of data, like base functions and associated
484  * methods to access those data on rows
485  * @param col_data structure of data, like base functions and associated
486  * methods to access those data on rows
487  * @return error code
488  */
489  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
490  EntityType col_type,
491  EntitiesFieldData::EntData &row_data,
492  EntitiesFieldData::EntData &col_data) {
493  MoFEMFunctionBegin;
494 
495  if (row_type == MBVERTEX)
496  MoFEMFunctionReturnHot(0);
497  if (col_type == MBVERTEX)
498  MoFEMFunctionReturnHot(0);
499 
500  const int nb_row_dofs = row_data.getN().size2() / 3;
501  if (nb_row_dofs == 0)
502  MoFEMFunctionReturnHot(0);
503  const int nb_col_dofs = col_data.getN().size2() / 3;
504  if (nb_col_dofs == 0)
505  MoFEMFunctionReturnHot(0);
506  entityLocalMatrix.resize(nb_row_dofs, nb_col_dofs, false);
507  entityLocalMatrix.clear();
508 
509  if (nb_row_dofs != static_cast<int>(row_data.getFieldData().size())) {
510  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
511  "Number of base functions and DOFs on entity is different on "
512  "rows %d!=%d",
513  nb_row_dofs, row_data.getFieldData().size());
514  }
515  if (nb_col_dofs != static_cast<int>(col_data.getFieldData().size())) {
516  SETERRQ2(
517  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
518  "Number of base functions and DOFs on entity is different on cols",
519  nb_col_dofs, col_data.getFieldData().size());
520  }
521 
522  MatrixDouble row_curl_mat, col_curl_mat;
523  FTensor::Index<'i', 3> i;
524  FTensor::Index<'j', 3> j;
525  FTensor::Index<'k', 3> k;
526 
527  const double c0 = 1. / blockData.mU;
528  const int nb_gauss_pts = row_data.getN().size1();
529  auto t_row_curl_base = row_data.getFTensor2DiffN<3, 3>();
530 
531  for (int gg = 0; gg != nb_gauss_pts; gg++) {
532 
533  // get integration weight scaled by volume
534  double w = getGaussPts()(3, gg) * getVolume();
535 
536  FTensor::Tensor1<double, 3> t_row_curl;
537  for (int aa = 0; aa != nb_row_dofs; aa++) {
538  t_row_curl(i) = levi_civita(j, i, k) * t_row_curl_base(j, k);
539 
540  FTensor::Tensor0<double *> t_local_mat(&entityLocalMatrix(aa, 0), 1);
541  FTensor::Tensor1<double, 3> t_col_curl;
542 
543  auto t_col_curl_base = col_data.getFTensor2DiffN<3, 3>(gg, 0);
544  for (int bb = 0; bb != nb_col_dofs; bb++) {
545  t_col_curl(i) = levi_civita(j, i, k) * t_col_curl_base(j, k);
546  t_local_mat += c0 * w * t_row_curl(i) * t_col_curl(i);
547  ++t_local_mat;
548  ++t_col_curl_base;
549  }
550 
551  ++t_row_curl_base;
552  }
553  }
554 
555  CHKERR MatSetValues(blockData.A, nb_row_dofs, &row_data.getIndices()[0],
556  nb_col_dofs, &col_data.getIndices()[0],
557  &entityLocalMatrix(0, 0), ADD_VALUES);
558 
559  if (row_side != col_side || row_type != col_type) {
560  entityLocalMatrix = trans(entityLocalMatrix);
561  CHKERR MatSetValues(blockData.A, nb_col_dofs, &col_data.getIndices()[0],
562  nb_row_dofs, &row_data.getIndices()[0],
563  &entityLocalMatrix(0, 0), ADD_VALUES);
564  }
565 
566  MoFEMFunctionReturn(0);
567  }
568  };
569 
570  /** \brief calculate and assemble stabilization matrix
571  * \ingroup maxwell_element
572 
573  \f[
574  \mathbf{A} = \int_\Omega \epsilon \mathbf{u} \cdot \mathbf{v}
575  \textrm{d}\Omega \f] where \f$\epsilon\f$ is regularization parameter.
576 
577  */
578  struct OpStab
579  : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
580 
581  BlockData &blockData;
582  OpStab(BlockData &data)
583  : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
584  data.fieldName, UserDataOperator::OPROWCOL),
585  blockData(data) {
586  sYmm = true;
587  }
588 
589  MatrixDouble entityLocalMatrix;
590 
591  /**
592  * \brief integrate matrix
593  * @param row_side local number of entity on element for row of the matrix
594  * @param col_side local number of entity on element for col of the matrix
595  * @param row_type type of row entity (EDGE/TRIANGLE/TETRAHEDRON)
596  * @param col_type type of col entity (EDGE/TRIANGLE/TETRAHEDRON)
597  * @param row_data structure of data, like base functions and associated
598  * methods to access those data on rows
599  * @param col_data structure of data, like base functions and associated
600  * methods to access those data on rows
601  * @return error code
602  */
603  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
604  EntityType col_type,
605  EntitiesFieldData::EntData &row_data,
606  EntitiesFieldData::EntData &col_data) {
607  MoFEMFunctionBegin;
608 
609  if (row_type == MBVERTEX)
610  MoFEMFunctionReturnHot(0);
611  if (col_type == MBVERTEX)
612  MoFEMFunctionReturnHot(0);
613 
614  const int nb_row_dofs = row_data.getN().size2() / 3;
615  if (nb_row_dofs == 0)
616  MoFEMFunctionReturnHot(0);
617  const int nb_col_dofs = col_data.getN().size2() / 3;
618  if (nb_col_dofs == 0)
619  MoFEMFunctionReturnHot(0);
620  entityLocalMatrix.resize(nb_row_dofs, nb_col_dofs, false);
621  entityLocalMatrix.clear();
622 
623  if (nb_row_dofs != static_cast<int>(row_data.getFieldData().size())) {
624  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
625  "Number of base functions and DOFs on entity is different on "
626  "rows %d!=%d",
627  nb_row_dofs, row_data.getFieldData().size());
628  }
629  if (nb_col_dofs != static_cast<int>(col_data.getFieldData().size())) {
630  SETERRQ2(
631  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
632  "Number of base functions and DOFs on entity is different on cols",
633  nb_col_dofs, col_data.getFieldData().size());
634  }
635 
636  MatrixDouble row_curl_mat, col_curl_mat;
637  FTensor::Index<'i', 3> i;
638 
639  const double c0 = 1. / blockData.mU;
640  const double c1 = blockData.ePsilon * c0;
641  const int nb_gauss_pts = row_data.getN().size1();
642 
643  for (int gg = 0; gg != nb_gauss_pts; gg++) {
644 
645  // get integration weight scaled by volume
646  double w = getGaussPts()(3, gg) * getVolume();
647 
648  FTensor::Tensor1<const double *, 3> t_row_base(
649  &row_data.getVectorN<3>(gg)(0, HVEC0),
650  &row_data.getVectorN<3>(gg)(0, HVEC1),
651  &row_data.getVectorN<3>(gg)(0, HVEC2), 3);
652 
653  for (int aa = 0; aa != nb_row_dofs; aa++) {
654  FTensor::Tensor0<double *> t_local_mat(&entityLocalMatrix(aa, 0), 1);
655  FTensor::Tensor1<const double *, 3> t_col_base(
656  &col_data.getVectorN<3>(gg)(0, HVEC0),
657  &col_data.getVectorN<3>(gg)(0, HVEC1),
658  &col_data.getVectorN<3>(gg)(0, HVEC2), 3);
659  for (int bb = 0; bb != nb_col_dofs; bb++) {
660  t_local_mat += c1 * w * t_row_base(i) * t_col_base(i);
661  ++t_col_base;
662  ++t_local_mat;
663  }
664  ++t_row_base;
665  }
666  }
667 
668  // cerr << entityLocalMatrix << endl;
669  // cerr << endl;
670 
671  CHKERR MatSetValues(blockData.A, nb_row_dofs, &row_data.getIndices()[0],
672  nb_col_dofs, &col_data.getIndices()[0],
673  &entityLocalMatrix(0, 0), ADD_VALUES);
674 
675  if (row_side != col_side || row_type != col_type) {
676  entityLocalMatrix = trans(entityLocalMatrix);
677  CHKERR MatSetValues(blockData.A, nb_col_dofs, &col_data.getIndices()[0],
678  nb_row_dofs, &row_data.getIndices()[0],
679  &entityLocalMatrix(0, 0), ADD_VALUES);
680  }
681 
682  MoFEMFunctionReturn(0);
683  }
684  };
685 
686  /** \brief calculate essential boundary conditions
687  * \ingroup maxwell_element
688 
689  \f[
690  \mathbf{F} = \int_{\partial\Omega} \ \mathbf{u} \cdot \mathbf{j}_i
691  \textrm{d}{\partial\Omega} \f] where \f$\mathbf{j}_i\f$ is current density
692  function.
693 
694  Here simple current on coil is hard coded. In future more general
695  implementation is needed.
696 
697  */
698  struct OpNaturalBC
699  : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
700 
701  BlockData &blockData;
702 
703  OpNaturalBC(BlockData &data)
704  : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
705  data.fieldName, UserDataOperator::OPROW),
706  blockData(data) {}
707 
708  VectorDouble naturalBC;
709 
710  /**
711  * \brief integrate matrix
712  * \ingroup maxwell_element
713  * @param row_side local number of entity on element for row of the matrix
714  * @param row_type type of row entity (EDGE/TRIANGLE/TETRAHEDRON)
715  * @param row_data structure of data, like base functions and associated
716  * methods to access those data on rows
717  * @return error code
718  */
719  MoFEMErrorCode doWork(int row_side, EntityType row_type,
720  EntitiesFieldData::EntData &row_data) {
721  MoFEMFunctionBegin;
722 
723  if (row_type == MBVERTEX)
724  MoFEMFunctionReturnHot(0);
725 
726  const int nb_row_dofs = row_data.getN().size2() / 3;
727  if (nb_row_dofs == 0)
728  MoFEMFunctionReturnHot(0);
729  naturalBC.resize(nb_row_dofs, false);
730  naturalBC.clear();
731 
732  FTensor::Index<'i', 3> i;
733 
734  const int nb_gauss_pts = row_data.getN().size1();
735  auto t_row_base = row_data.getFTensor1N<3>();
736 
737  for (int gg = 0; gg != nb_gauss_pts; gg++) {
738 
739  // get integration weight scaled by volume
740  double area;
741  area = norm_2(getNormalsAtGaussPts(gg)) * 0.5;
742  double w = getGaussPts()(2, gg) * area;
743 
744  // Current is on surface where natural bc are applied. It is set that
745  // current is in XY plane, circular, around the coil.
746  const double x = getCoordsAtGaussPts()(gg, 0);
747  const double y = getCoordsAtGaussPts()(gg, 1);
748  const double r = sqrt(x * x + y * y);
749  FTensor::Tensor1<double, 3> t_j;
750  t_j(0) = -y / r;
751  t_j(1) = +x / r;
752  t_j(2) = 0;
753 
754  // double a = t_j(i) * t_tangent1(i);
755  // double b = t_j(i) * t_tangent2(i);
756  // t_j(i) = a * t_tangent1(i) + b * t_tangent2(i);
757 
758  // ++t_tangent1;
759  // ++t_tangent2;
760 
761  FTensor::Tensor0<double *> t_f(&naturalBC[0]);
762  for (int aa = 0; aa != nb_row_dofs; aa++) {
763  t_f += w * t_row_base(i) * t_j(i);
764  ++t_row_base;
765  ++t_f;
766  }
767  }
768 
769  CHKERR VecSetValues(blockData.F, row_data.getIndices().size(),
770  &row_data.getIndices()[0], &naturalBC[0], ADD_VALUES);
771 
772  MoFEMFunctionReturn(0);
773  }
774  };
775 
776  /** \brief calculate and assemble CurlCurl matrix
777  * \ingroup maxwell_element
778  */
779  struct OpPostProcessCurl
780  : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
781 
782  BlockData &blockData;
783  moab::Interface &postProcMesh;
784  std::vector<EntityHandle> &mapGaussPts;
785 
786  OpPostProcessCurl(BlockData &data, moab::Interface &post_proc_mesh,
787  std::vector<EntityHandle> &map_gauss_pts)
788  : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
789  data.fieldName, UserDataOperator::OPROW),
790  blockData(data), postProcMesh(post_proc_mesh),
791  mapGaussPts(map_gauss_pts) {}
792 
793  MoFEMErrorCode doWork(int row_side, EntityType row_type,
794  EntitiesFieldData::EntData &row_data) {
795  MoFEMFunctionBegin;
796 
797  if (row_type == MBVERTEX)
798  MoFEMFunctionReturnHot(0);
799 
800  Tag th;
801  double def_val[] = {0, 0, 0};
802  CHKERR postProcMesh.tag_get_handle("MAGNETIC_INDUCTION_FIELD", 3,
803  MB_TYPE_DOUBLE, th,
804  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
805  const int nb_row_dofs = row_data.getN().size2() / 3;
806  if (nb_row_dofs == 0)
807  MoFEMFunctionReturnHot(0);
808  const void *tags_ptr[mapGaussPts.size()];
809 
810  if (nb_row_dofs != row_data.getFieldData().size())
811  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
812  "Wrong number of base functions and DOFs %d != %d",
813  nb_row_dofs, row_data.getFieldData().size());
814 
815  FTensor::Index<'i', 3> i;
816  FTensor::Index<'j', 3> j;
817  FTensor::Index<'k', 3> k;
818  const int nb_gauss_pts = row_data.getN().size1();
819  if (nb_gauss_pts != static_cast<int>(mapGaussPts.size())) {
820  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
821  "Inconsistency number of dofs %d!=%d", nb_gauss_pts,
822  mapGaussPts.size());
823  }
824  CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0],
825  mapGaussPts.size(), tags_ptr);
826  auto t_curl_base = row_data.getFTensor2DiffN<3, 3>();
827  for (int gg = 0; gg != nb_gauss_pts; gg++) {
828  // get pointer to tag values on entity (i.e. vertex on refined
829  // post-processing mesh)
830  double *ptr = &((double *)tags_ptr[gg])[0];
831  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_curl(ptr, &ptr[1],
832  &ptr[2]);
833  // calculate curl value
834  auto t_dof = row_data.getFTensor0FieldData();
835  for (int aa = 0; aa != nb_row_dofs; aa++) {
836  t_curl(i) += t_dof * (levi_civita(j, i, k) * t_curl_base(j, k));
837  ++t_curl_base;
838  ++t_dof;
839  }
840  }
841  MoFEMFunctionReturn(0);
842  }
843  };
844 };
845 
846 #endif //__MAGNETICELEMENT_HPP__
847 
848 /******************************************************************************
849  * \defgroup maxwell_element Magnetic/Maxwell element
850  * \ingroup user_modules
851  ******************************************************************************/
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MagneticElement::OpNaturalBC::naturalBC
VectorDouble naturalBC
Definition: MagneticElement.hpp:708
MoFEM::VolumeElementForcesAndSourcesCore::VolumeElementForcesAndSourcesCore
VolumeElementForcesAndSourcesCore(Interface &m_field, const EntityType type=MBTET)
Definition: VolumeElementForcesAndSourcesCore.cpp:9
MagneticElement::BlockData
data structure storing material constants, model parameters, matrices, etc.
Definition: MagneticElement.hpp:64
MagneticElement::BlockData::~BlockData
~BlockData()
Definition: MagneticElement.hpp:91
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
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.
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1644
FTensor::Tensor1< double, 3 >
EntityHandle
MagneticElement::BlockData::D
Vec D
Definition: MagneticElement.hpp:86
MoFEM::addHOOpsVol
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Definition: HODataOperators.hpp:674
MoFEM::FaceElementForcesAndSourcesCore::meshPositionsFieldName
std::string meshPositionsFieldName
Definition: FaceElementForcesAndSourcesCore.hpp:29
MagneticElement::BlockData::mU
double mU
magnetic constant N / A2
Definition: MagneticElement.hpp:72
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
MagneticElement::BlockData::dM
DM dM
Definition: MagneticElement.hpp:84
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MagneticElement::destroyProblem
MoFEMErrorCode destroyProblem()
destroy Distributed mesh manager
Definition: MagneticElement.hpp:307
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
MagneticElement::OpStab::entityLocalMatrix
MatrixDouble entityLocalMatrix
Definition: MagneticElement.hpp:589
MoFEM::DMoFEMMeshToLocalVector
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
MagneticElement::BlockData::essentialBc
Range essentialBc
Definition: MagneticElement.hpp:79
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition: Levi_Civita.hpp:617
MagneticElement::MagneticElement
MagneticElement(MoFEM::Interface &m_field)
Definition: MagneticElement.hpp:56
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
HVEC1
@ HVEC1
Definition: definitions.h:199
sdf.r
int r
Definition: sdf.py:8
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::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MagneticElement::OpPostProcessCurl::blockData
BlockData & blockData
Definition: MagneticElement.hpp:782
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MagneticElement::BlockData::fieldName
const string fieldName
Definition: MagneticElement.hpp:67
MagneticElement::BlockData::BlockData
BlockData()
Definition: MagneticElement.hpp:88
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
MagneticElement::BlockData::F
Vec F
Definition: MagneticElement.hpp:86
MagneticElement::OpCurlCurl::OpCurlCurl
OpCurlCurl(BlockData &data)
Definition: MagneticElement.hpp:468
MagneticElement::OpNaturalBC::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
integrate matrix
Definition: MagneticElement.hpp:719
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
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::get_moab
virtual moab::Interface & get_moab()=0
MagneticElement::getEssentialBc
MoFEMErrorCode getEssentialBc()
get essential boundary conditions (only homogenous case is considered)
Definition: MagneticElement.hpp:121
MagneticElement::OpStab
calculate and assemble stabilization matrix
Definition: MagneticElement.hpp:578
MagneticElement::createFields
MoFEMErrorCode createFields()
build problem data structures
Definition: MagneticElement.hpp:158
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
MagneticElement::VolumeFE::VolumeFE
VolumeFE(MoFEM::Interface &m_field)
Definition: MagneticElement.hpp:34
MagneticElement::BlockData::naturalBc
Range naturalBc
Definition: MagneticElement.hpp:76
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:489
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
MagneticElement::OpPostProcessCurl::OpPostProcessCurl
OpPostProcessCurl(BlockData &data, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts)
Definition: MagneticElement.hpp:786
MagneticElement::TriFE::TriFE
TriFE(MoFEM::Interface &m_field)
Definition: MagneticElement.hpp:51
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MagneticElement::BlockData::oRder
int oRder
approximation order
Definition: MagneticElement.hpp:81
MagneticElement::VolumeFE
definition of volume element
Definition: MagneticElement.hpp:33
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MagneticElement::OpPostProcessCurl::mapGaussPts
std::vector< EntityHandle > & mapGaussPts
Definition: MagneticElement.hpp:784
MagneticElement::BlockData::feNaturalBCName
const string feNaturalBCName
Definition: MagneticElement.hpp:69
MagneticElement::postProcessResults
MoFEMErrorCode postProcessResults()
post-process results, i.e. save solution on the mesh
Definition: MagneticElement.hpp:406
MoFEM::DMMoFEMCreateMoFEM
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMoFEM.cpp:114
MagneticElement::OpNaturalBC::blockData
BlockData & blockData
Definition: MagneticElement.hpp:701
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::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MagneticElement::OpStab::OpStab
OpStab(BlockData &data)
Definition: MagneticElement.hpp:582
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:108
FaceElementForcesAndSourcesCore
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
FTensor::Index< 'i', 3 >
MagneticElement::TriFE::getRule
int getRule(int order)
Definition: MagneticElement.hpp:53
MagneticElement::blockData
BlockData blockData
Definition: MagneticElement.hpp:94
MagneticElement::OpPostProcessCurl::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Definition: MagneticElement.hpp:793
MagneticElement::OpNaturalBC::OpNaturalBC
OpNaturalBC(BlockData &data)
Definition: MagneticElement.hpp:703
MagneticElement::createProblem
MoFEMErrorCode createProblem()
create problem
Definition: MagneticElement.hpp:278
Range
MagneticElement::TriFE
define surface element
Definition: MagneticElement.hpp:50
MagneticElement::mField
MoFEM::Interface & mField
Definition: MagneticElement.hpp:30
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
FTensor::Tensor0
Definition: Tensor0.hpp:16
_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:71
MagneticElement::OpPostProcessCurl
calculate and assemble CurlCurl matrix
Definition: MagneticElement.hpp:779
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MagneticElement::BlockData::ePsilon
double ePsilon
regularization paramater
Definition: MagneticElement.hpp:73
MagneticElement::OpCurlCurl::blockData
BlockData & blockData
Definition: MagneticElement.hpp:467
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MagneticElement::getNaturalBc
MoFEMErrorCode getNaturalBc()
get natural boundary conditions
Definition: MagneticElement.hpp:101
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::DMMoFEMGetProblemPtr
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:426
MagneticElement::BlockData::feName
const string feName
Definition: MagneticElement.hpp:68
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
sdf_wavy_2d.w
int w
Definition: sdf_wavy_2d.py:6
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:603
MagneticElement::OpNaturalBC
calculate essential boundary conditions
Definition: MagneticElement.hpp:698
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MagneticElement::BlockData::A
Mat A
Definition: MagneticElement.hpp:85
MagneticElement::createElements
MoFEMErrorCode createElements()
create finite elements
Definition: MagneticElement.hpp:234
MagneticElement
Implementation of magneto-static problem (basic Implementation)
Definition: MagneticElement.hpp:28
MagneticElement::OpStab::blockData
BlockData & blockData
Definition: MagneticElement.hpp:581
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
HVEC2
@ HVEC2
Definition: definitions.h:199
MagneticElement::OpPostProcessCurl::postProcMesh
moab::Interface & postProcMesh
Definition: MagneticElement.hpp:783
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::Problem::getName
auto getName() const
Definition: ProblemsMultiIndices.hpp:372
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
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:586
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MagneticElement::OpCurlCurl
calculate and assemble CurlCurl matrix
Definition: MagneticElement.hpp:464
MagneticElement::VolumeFE::getRule
int getRule(int order)
Definition: MagneticElement.hpp:36
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.
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
tol
double tol
Definition: mesh_smoothing.cpp:26
HVEC0
@ HVEC0
Definition: definitions.h:199
MagneticElement::solveProblem
MoFEMErrorCode solveProblem(const bool regression_test=false)
solve problem
Definition: MagneticElement.hpp:319
MagneticElement::~MagneticElement
virtual ~MagneticElement()=default
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1123
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
MagneticElement::OpCurlCurl::entityLocalMatrix
MatrixDouble entityLocalMatrix
Definition: MagneticElement.hpp:475
Generated by Doxygen 1.8.17 and hosted at University of Glasgow