v0.14.0
continuity_check_on_skeleton_with_simple_2d_for_h1.cpp
Go to the documentation of this file.
1 /**
2  * \file continuity_check_on_skeleton_with_simple_2d_for_h1.cpp
3  * \ingroup mofem_simple_interface
4  \example continuity_check_on_skeleton_with_simple_2d_for_h1.cpp
5  *
6  * \brief Integration on skeleton for 2d
7  *
8  * Teting integration on skeleton and checking of continuity of hcurl space on
9  * edges.
10  */
11 
12 #include <MoFEM.hpp>
13 using namespace MoFEM;
14 
15 static char help[] = "...\n\n";
16 
21 
22 struct CommonData {
26 };
27 
28 struct SkeletonFE : public EdgeEleOp {
29 
30  struct OpFaceSide : public FaceEleOnSideOp {
31 
33  OpFaceSide(CommonData &elem_data)
34  : FaceEleOnSideOp("FIELD", UserDataOperator::OPROW),
35  elemData(elem_data) {}
36 
37  MoFEMErrorCode doWork(int side, EntityType type,
40 
41  if (type == MBEDGE && side == getEdgeSideNumber()) {
42 
43  MatrixDouble diff = getCoordsAtGaussPts() - getEdgeCoordsAtGaussPts();
44 
45  const double eps = 1e-12;
46  if (norm_inf(diff) > eps) {
47  MOFEM_LOG_ATTRIBUTES("ATOM",
50  MOFEM_LOG("ATOM", Sev::error)
51  << "Quad coords: " << getCoordsAtGaussPts();
52  MOFEM_LOG("ATOM", Sev::error)
53  << "Edge coords: " << getEdgeCoordsAtGaussPts();
54  MOFEM_LOG("ATOM", Sev::error) << "Difference: " << diff;
55  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
56  "Coordinates at integration pts are different");
57  }
58 
59  const size_t nb_dofs = data.getN().size2();
60  const size_t nb_integration_pts = data.getN().size1();
61 
62  auto t_base = data.getFTensor0N();
63  MatrixDouble *ptr_dot_elem_data = nullptr;
64  if (getFEMethod()->nInTheLoop == 0)
65  ptr_dot_elem_data = &elemData.dotEleLeft;
66  else
67  ptr_dot_elem_data = &elemData.dotEleRight;
68  MatrixDouble &dot_elem_data = *ptr_dot_elem_data;
69  dot_elem_data.resize(nb_integration_pts, nb_dofs, false);
70 
71  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
72  for (size_t bb = 0; bb != nb_dofs; ++bb) {
73  dot_elem_data(gg, bb) = t_base;
74  ++t_base;
75  }
76  }
77  }
79  }
80  };
81 
84 
85  SkeletonFE(MoFEM::Interface &m_field, CommonData &elem_data)
86  : EdgeEle::UserDataOperator("FIELD", UserDataOperator::OPROW),
87  faceSideFe(m_field), elemData(elem_data) {
88  faceSideFe.getOpPtrVector().push_back(new SkeletonFE::OpFaceSide(elemData));
89  }
90 
91  MoFEMErrorCode doWork(int side, EntityType type,
93 
95  if (type == MBEDGE) {
96 
97  const size_t nb_dofs = data.getN().size2();
98  const size_t nb_integration_pts = data.getN().size1();
99 
100  MOFEM_LOG("ATOM", Sev::noisy)
101  << "Cords at integration points" << getCoordsAtGaussPts();
102 
103  auto t_base = data.getFTensor0N();
104  elemData.dotEdge.resize(nb_integration_pts, nb_dofs, false);
105  elemData.dotEleLeft.resize(nb_integration_pts, 0, false);
106  elemData.dotEleRight.resize(nb_integration_pts, 0, false);
107 
108  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
109  for (size_t bb = 0; bb != nb_dofs; ++bb) {
110  elemData.dotEdge(gg, bb) = t_base;
111  ++t_base;
112  }
113  }
114 
115  CHKERR loopSideFaces("dFE", faceSideFe);
116 
117  auto check_continuity_of_base = [&](auto &vol_dot_data) {
119 
120  if (vol_dot_data.size1() != elemData.dotEdge.size1())
121  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
122  "Inconsistent number of integration points");
123 
124  if (vol_dot_data.size2() != elemData.dotEdge.size2())
125  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
126  "Inconsistent number of base functions");
127  const double eps = 1e-12;
128  for (size_t gg = 0; gg != vol_dot_data.size1(); ++gg)
129  for (size_t bb = 0; bb != vol_dot_data.size2(); ++bb) {
130  const double error =
131  std::abs(vol_dot_data(gg, bb) - elemData.dotEdge(gg, bb));
132  if (error > eps)
133  SETERRQ4(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
134  "Inconsistency (%d, %d) %3.4e != %3.4e", gg, bb,
135  vol_dot_data(gg, bb), elemData.dotEdge(gg, bb));
136  else
137  MOFEM_LOG("ATOM", Sev::noisy) << "Ok";
138  }
140  };
141 
142  if (elemData.dotEleLeft.size2() != 0)
143  CHKERR check_continuity_of_base(elemData.dotEleLeft);
144  else if (elemData.dotEleRight.size2() != 0)
145  CHKERR check_continuity_of_base(elemData.dotEleRight);
146  }
148  }
149 };
150 
151 int main(int argc, char *argv[]) {
152 
153  // initialize petsc
154  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
155 
156  try {
157 
158  // Create MoAB database
159  moab::Core moab_core;
160  moab::Interface &moab = moab_core;
161 
162  // Create MoFEM database and link it to MoAB
163  MoFEM::Core mofem_core(moab);
164  MoFEM::Interface &m_field = mofem_core;
165 
166  auto core_log = logging::core::get();
167  core_log->add_sink(
169  LogManager::setLog("ATOM");
170 
171  // Register DM Manager
172  DMType dm_name = "DMMOFEM";
173  CHKERR DMRegister_MoFEM(dm_name);
174 
175  // Simple interface
176  Simple *simple_interface;
177  CHKERR m_field.getInterface(simple_interface);
178  {
179  // get options from command line
180  CHKERR simple_interface->getOptions();
181  // Add ghost cells needed for evaluation on skeleton
182  simple_interface->getAddSkeletonFE() = true;
183 
184  // load mesh file
185  CHKERR simple_interface->loadFile("");
186 
187  auto get_base = []() -> FieldApproximationBase {
188  enum bases { AINSWORTH, DEMKOWICZ, LASTBASEOP };
189  const char *list_bases[] = {"ainsworth", "demkowicz"};
190  PetscBool flg;
191  PetscInt choice_base_value = AINSWORTH;
192  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
193  LASTBASEOP, &choice_base_value, &flg);
194  if (flg == PETSC_TRUE) {
196  if (choice_base_value == AINSWORTH)
198  else if (choice_base_value == DEMKOWICZ)
199  base = DEMKOWICZ_JACOBI_BASE;
200  return base;
201  }
202  return LASTBASE;
203  };
204 
205  // add fields
206  auto base = get_base();
207  CHKERR simple_interface->addDomainField("FIELD", H1, base, 1);
208  // CHKERR simple_interface->addDomainField("TEST_FIELD", L2,
209  // AINSWORTH_LEGENDRE_BASE, 1);
210  CHKERR simple_interface->addSkeletonField("FIELD", H1, base, 1);
211  // set fields order
212  CHKERR simple_interface->setFieldOrder("FIELD", 2);
213  // CHKERR simple_interface->setFieldOrder("TEST_FIELD", 1);
214  // setup problem
215  CHKERR simple_interface->setUp();
216  // get dm
217  auto dm = simple_interface->getDM();
218 
219  // create elements
220  CommonData elem_data;
221  boost::shared_ptr<EdgeEle> skeleton_fe =
222  boost::shared_ptr<EdgeEle>(new EdgeEle(m_field));
223  skeleton_fe->getOpPtrVector().push_back(
224  new SkeletonFE(m_field, elem_data));
225 
226  // iterate skeleton finite elements
227  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getSkeletonFEName(),
228  skeleton_fe);
229  }
230  }
231  CATCH_ERRORS;
232 
233  // finish work cleaning memory, getting statistics, etc.
235 
236  return 0;
237 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::EdgeElementForcesAndSourcesCore
Edge finite element.
Definition: EdgeElementForcesAndSourcesCore.hpp:30
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
CommonData::dotEdge
MatrixDouble dotEdge
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:23
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
LASTBASE
@ LASTBASE
Definition: definitions.h:69
SkeletonFE::OpFaceSide
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:30
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
CommonData::dotEleLeft
MatrixDouble dotEleLeft
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:24
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
SkeletonFE::faceSideFe
FaceEleOnSide faceSideFe
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:83
MoFEM::Simple::getSkeletonFEName
const std::string getSkeletonFEName() const
Get the Skeleton FE Name.
Definition: Simple.hpp:355
MOFEM_LOG_ATTRIBUTES
#define MOFEM_LOG_ATTRIBUTES(channel, bit)
Add attributes to channel.
Definition: LogManager.hpp:296
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
MoFEM.hpp
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
SkeletonFE::OpFaceSide::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:37
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
SkeletonFE
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:28
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
SkeletonFE::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:91
MoFEM::Simple::addSkeletonField
MoFEMErrorCode addSkeletonField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on skeleton.
Definition: Simple.cpp:300
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1489
MOFEM_LOG_FUNCTION
#define MOFEM_LOG_FUNCTION()
Set scope.
Definition: LogManager.hpp:325
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEM::Simple::getAddSkeletonFE
bool & getAddSkeletonFE()
Get the addSkeletonFE.
Definition: Simple.hpp:425
MoFEM::LogManager::BitLineID
@ BitLineID
Definition: LogManager.hpp:48
MoFEM::LogManager::BitScope
@ BitScope
Definition: LogManager.hpp:49
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:667
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: FaceElementForcesAndSourcesCore.hpp:86
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
EdgeEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:18
convert.type
type
Definition: convert.py:64
CommonData::dotEleRight
MatrixDouble dotEleRight
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:25
MoFEM::Simple::addDomainField
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
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
SkeletonFE::OpFaceSide::elemData
CommonData & elemData
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:32
SkeletonFE::OpFaceSide::OpFaceSide
OpFaceSide(CommonData &elem_data)
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:33
MoFEM::LogManager::getStrmSelf
static boost::shared_ptr< std::ostream > getStrmSelf()
Get the strm self object.
Definition: LogManager.cpp:340
MoFEM::ForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: ForcesAndSourcesCore.hpp:482
help
static char help[]
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:15
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
MoFEM::EdgeElementForcesAndSourcesCore::UserDataOperator
default operator for EDGE element
Definition: EdgeElementForcesAndSourcesCore.hpp:68
SkeletonFE::SkeletonFE
SkeletonFE(MoFEM::Interface &m_field, CommonData &elem_data)
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:85
MoFEM::CoreTmp< 0 >::Initialize
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
main
int main(int argc, char *argv[])
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:151
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: EntitiesFieldData.hpp:1305
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1102
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
CommonData
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:22
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:203
MoFEM::FaceElementForcesAndSourcesCoreOnSide::UserDataOperator
default operator for Face element
Definition: FaceElementForcesAndSourcesCoreOnSide.hpp:92
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
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:416
MoFEM::FaceElementForcesAndSourcesCoreOnSide
Base face element used to integrate on skeleton.
Definition: FaceElementForcesAndSourcesCoreOnSide.hpp:19
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
SkeletonFE::elemData
CommonData & elemData
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:82