v0.10.0
dynamic_elastic.cpp
Go to the documentation of this file.
1 /**
2  * \file dynamic_elastic.cpp
3  * \example dynamic_elastic.cpp
4  *
5  * Plane stress elastic dynamic problem
6  *
7  */
8 
9 /* This file is part of MoFEM.
10  * MoFEM is free software: you can redistribute it and/or modify it under
11  * the terms of the GNU Lesser General Public License as published by the
12  * Free Software Foundation, either version 3 of the License, or (at your
13  * option) any later version.
14  *
15  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
16  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18  * License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
22 
23 #include <MoFEM.hpp>
24 
25 using namespace MoFEM;
26 
27 template <int DIM> struct ElementsAndOps {};
28 
29 template <> struct ElementsAndOps<2> {
33 };
34 
35 template <> struct ElementsAndOps<3> {
39 };
40 
41 constexpr int SPACE_DIM = 2; //< Space dimension of problem, mesh
42 
47 
49  GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
55  GAUSS>::OpSource<1, SPACE_DIM>;
58 
59 constexpr double rho = 1;
60 constexpr double omega = 2.4;
61 constexpr double young_modulus = 1;
62 constexpr double poisson_ratio = 0.25;
63 constexpr double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
64 constexpr double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
65 
66 #include <OpPostProcElastic.hpp>
67 using namespace Tutorial;
68 
69 static double *ts_time_ptr;
70 static double *ts_aa_ptr;
71 
72 struct Example {
73 
74  Example(MoFEM::Interface &m_field) : mField(m_field) {}
75 
76  MoFEMErrorCode runProblem();
77 
78 private:
79  MoFEM::Interface &mField;
80 
81  MoFEMErrorCode readMesh();
82  MoFEMErrorCode setupProblem();
83  MoFEMErrorCode createCommonData();
84  MoFEMErrorCode boundaryCondition();
85  MoFEMErrorCode assembleSystem();
86  MoFEMErrorCode solveSystem();
87  MoFEMErrorCode outputResults();
88  MoFEMErrorCode checkResults();
89 
91  boost::shared_ptr<MatrixDouble> matGradPtr;
92  boost::shared_ptr<MatrixDouble> matStrainPtr;
93  boost::shared_ptr<MatrixDouble> matStressPtr;
94  boost::shared_ptr<MatrixDouble> matAccelerationPtr;
95  boost::shared_ptr<MatrixDouble> matInertiaPtr;
96  boost::shared_ptr<MatrixDouble> matDPtr;
97 };
98 
99 //! [Create common data]
102 
103  auto set_matrial_stiffens = [&]() {
108  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
110  constexpr double A =
111  (SPACE_DIM == 2) ? 2 * shear_modulus_G /
112  (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
113  : 1;
114  auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*matDPtr);
115  t_D(i, j, k, l) = 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
116  A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
117  t_kd(i, j) * t_kd(k, l);
119  };
120 
121  matGradPtr = boost::make_shared<MatrixDouble>();
122  matStrainPtr = boost::make_shared<MatrixDouble>();
123  matStressPtr = boost::make_shared<MatrixDouble>();
124  matAccelerationPtr = boost::make_shared<MatrixDouble>();
125  matInertiaPtr = boost::make_shared<MatrixDouble>();
126  matDPtr = boost::make_shared<MatrixDouble>();
127 
128  constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
129  matDPtr->resize(size_symm * size_symm, 1);
130 
131  CHKERR set_matrial_stiffens();
132 
134 }
135 //! [Create common data]
136 
137 //! [Run problem]
140  CHKERR readMesh();
141  CHKERR setupProblem();
142  CHKERR createCommonData();
143  CHKERR boundaryCondition();
144  CHKERR assembleSystem();
145  CHKERR solveSystem();
146  CHKERR outputResults();
147  CHKERR checkResults();
149 }
150 //! [Run problem]
151 
152 //! [Read mesh]
155  auto simple = mField.getInterface<Simple>();
156  CHKERR simple->getOptions();
157  CHKERR simple->loadFile();
159 }
160 //! [Read mesh]
161 
162 //! [Set up problem]
165  Simple *simple = mField.getInterface<Simple>();
166  // Add field
167  CHKERR simple->addDomainField("U", H1, AINSWORTH_BERNSTEIN_BEZIER_BASE,
168  SPACE_DIM);
169  int order = 3;
170  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
171  CHKERR simple->setFieldOrder("U", order);
172  CHKERR simple->setUp();
174 }
175 //! [Set up problem]
176 
177 //! [Boundary condition]
180 
181  auto fix_disp = [&](const std::string blockset_name) {
182  Range fix_ents;
184  if (it->getName().compare(0, blockset_name.length(), blockset_name) ==
185  0) {
186  CHKERR mField.get_moab().get_entities_by_handle(it->meshset, fix_ents,
187  true);
188  }
189  }
190  return fix_ents;
191  };
192 
193  auto remove_ents = [&](const Range &&ents, const int lo, const int hi) {
194  auto prb_mng = mField.getInterface<ProblemsManager>();
195  auto simple = mField.getInterface<Simple>();
197  Range verts;
198  CHKERR mField.get_moab().get_connectivity(ents, verts, true);
199  verts.merge(ents);
200  if (SPACE_DIM == 3) {
201  Range adj;
202  CHKERR mField.get_moab().get_adjacencies(ents, 1, false, adj,
203  moab::Interface::UNION);
204  verts.merge(adj);
205  };
206 
207  CHKERR prb_mng->removeDofsOnEntities(simple->getProblemName(), "U", verts,
208  lo, hi);
210  };
211 
212  CHKERR remove_ents(fix_disp("FIX_X"), 0, 0);
213  CHKERR remove_ents(fix_disp("FIX_Y"), 1, 1);
214  CHKERR remove_ents(fix_disp("FIX_Z"), 2, 2);
215  CHKERR remove_ents(fix_disp("FIX_ALL"), 0, 3);
216 
218 }
219 //! [Boundary condition]
220 
221 //! [Push operators to pipeline]
224  auto *simple = mField.getInterface<Simple>();
225  auto *pipeline_mng = mField.getInterface<PipelineManager>();
226 
227  if (SPACE_DIM == 2) {
228  pipeline_mng->getOpDomainLhsPipeline().push_back(
230  pipeline_mng->getOpDomainLhsPipeline().push_back(
232  pipeline_mng->getOpDomainRhsPipeline().push_back(
234  pipeline_mng->getOpDomainRhsPipeline().push_back(
236  }
237 
238  // Get pointer to U_tt shift in domain element
239  auto get_rho = [this](const double, const double, const double) {
240  auto *pipeline_mng = mField.getInterface<PipelineManager>();
241  auto &fe_domain_lhs = pipeline_mng->getDomainLhsFE();
242  return rho * fe_domain_lhs->ts_aa;
243  };
244 
245  auto get_body_force = [this](const double, const double, const double) {
246  auto *pipeline_mng = mField.getInterface<PipelineManager>();
247  auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
248  auto t_source = FTensor::Tensor1<double, SPACE_DIM>{0.1, 1.};
250  const auto time = fe_domain_rhs->ts_t;
251  t_source(i) *= sin(time * omega * M_PI);
252  return t_source;
253  };
254 
255  pipeline_mng->getOpDomainLhsPipeline().push_back(new OpK("U", "U", matDPtr));
256  pipeline_mng->getOpDomainLhsPipeline().push_back(
257  new OpMass("U", "U", get_rho));
258 
259  pipeline_mng->getOpDomainRhsPipeline().push_back(
260  new OpBodyForce("U", get_body_force));
261  pipeline_mng->getOpDomainRhsPipeline().push_back(
263  matGradPtr));
264  pipeline_mng->getOpDomainRhsPipeline().push_back(
265  new OpSymmetrizeTensor<SPACE_DIM>("U", matGradPtr, matStrainPtr));
266  pipeline_mng->getOpDomainRhsPipeline().push_back(
268  "U", matStrainPtr, matStressPtr, matDPtr));
269  pipeline_mng->getOpDomainRhsPipeline().push_back(
270  new OpInternalForce("U", matStressPtr));
271  pipeline_mng->getOpDomainRhsPipeline().push_back(
273  matAccelerationPtr));
274  pipeline_mng->getOpDomainRhsPipeline().push_back(
275  new OpScaleMatrix("U", rho, matAccelerationPtr, matInertiaPtr));
276  pipeline_mng->getOpDomainRhsPipeline().push_back(
277  new OpInertiaForce("U", matInertiaPtr));
278 
279  auto integration_rule = [](int, int, int approx_order) {
280  return 2 * approx_order;
281  };
282  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
283  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
284 
286 }
287 //! [Push operators to pipeline]
288 
289 /**
290  * @brief Monitor solution
291  *
292  * This functions is called by TS solver at the end of each step. It is used
293  * to output results to the hard drive.
294  */
295 struct Monitor : public FEMethod {
296  Monitor(SmartPetscObj<DM> dm, boost::shared_ptr<PostProcEle> post_proc)
297  : dM(dm), postProc(post_proc){};
300  constexpr int save_every_nth_step = 1;
301  if (ts_step % save_every_nth_step == 0) {
302  CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProc);
303  CHKERR postProc->writeFile(
304  "out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
305  }
307  }
308 
309 private:
311  boost::shared_ptr<PostProcEle> postProc;
312 };
313 
314 //! [Solve]
317  auto *simple = mField.getInterface<Simple>();
318  auto *pipeline_mng = mField.getInterface<PipelineManager>();
319 
320  auto dm = simple->getDM();
321  auto ts = pipeline_mng->createTS2();
322 
323  // Setup postprocessing
324  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
325  post_proc_fe->generateReferenceElementMesh();
326  if (SPACE_DIM) {
327  post_proc_fe->getOpPtrVector().push_back(
329  post_proc_fe->getOpPtrVector().push_back(new OpSetInvJacH1ForFace(invJac));
330  }
331  post_proc_fe->getOpPtrVector().push_back(
333  matGradPtr));
334  post_proc_fe->getOpPtrVector().push_back(
335  new OpSymmetrizeTensor<SPACE_DIM>("U", matGradPtr, matStrainPtr));
336  post_proc_fe->getOpPtrVector().push_back(
338  "U", matStrainPtr, matStressPtr, matDPtr));
339  post_proc_fe->getOpPtrVector().push_back(new OpPostProcElastic<SPACE_DIM>(
340  "U", post_proc_fe->postProcMesh, post_proc_fe->mapGaussPts, matStrainPtr,
341  matStressPtr));
342  post_proc_fe->addFieldValuesPostProc("U");
343 
344  // Add monitor to time solver
345  boost::shared_ptr<FEMethod> null_fe;
346  auto monitor_ptr = boost::make_shared<Monitor>(dm, post_proc_fe);
347  CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
348  null_fe, monitor_ptr);
349 
350  double ftime = 1;
351  CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
352  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
353 
354  auto T = smartCreateDMVector(simple->getDM());
355  auto TT = smartVectorDuplicate(T);
356  CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
357  SCATTER_FORWARD);
358 
359  CHKERR TS2SetSolution(ts, T, TT);
360  CHKERR TSSetFromOptions(ts);
361 
362  CHKERR TSSolve(ts, NULL);
363  CHKERR TSGetTime(ts, &ftime);
364 
365  PetscInt steps, snesfails, rejects, nonlinits, linits;
366  CHKERR TSGetTimeStepNumber(ts, &steps);
367  CHKERR TSGetSNESFailures(ts, &snesfails);
368  CHKERR TSGetStepRejections(ts, &rejects);
369  CHKERR TSGetSNESIterations(ts, &nonlinits);
370  CHKERR TSGetKSPIterations(ts, &linits);
371  MOFEM_LOG_C("EXAMPLE", Sev::inform,
372  "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits "
373  "%D, linits %D\n",
374  steps, rejects, snesfails, ftime, nonlinits, linits);
375 
377 }
378 //! [Solve]
379 
380 //! [Postprocess results]
383  PetscBool test_flg = PETSC_FALSE;
384  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test_flg, PETSC_NULL);
385  if (test_flg) {
386  auto *simple = mField.getInterface<Simple>();
387  auto T = smartCreateDMVector(simple->getDM());
388  CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
389  SCATTER_FORWARD);
390  double nrm2;
391  CHKERR VecNorm(T, NORM_2, &nrm2);
392  MOFEM_LOG("EXAMPLE", Sev::inform) << "Regression norm " << nrm2;
393  constexpr double regression_value = 1.09572;
394  if (fabs(nrm2 - regression_value) > 1e-2)
395  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
396  "Regression test faileed; wrong norm value.");
397  }
399 }
400 //! [Postprocess results]
401 
402 //! [Check]
406 }
407 //! [Check]
408 
409 static char help[] = "...\n\n";
410 
411 int main(int argc, char *argv[]) {
412 
413  // Initialisation of MoFEM/PETSc and MOAB data structures
414  const char param_file[] = "param_file.petsc";
415  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
416 
417  // Add logging channel for example
418  auto core_log = logging::core::get();
419  core_log->add_sink(
421  LogManager::setLog("EXAMPLE");
422  MOFEM_LOG_TAG("EXAMPLE", "example");
423 
424  try {
425 
426  //! [Register MoFEM discrete manager in PETSc]
427  DMType dm_name = "DMMOFEM";
428  CHKERR DMRegister_MoFEM(dm_name);
429  //! [Register MoFEM discrete manager in PETSc
430 
431  //! [Create MoAB]
432  moab::Core mb_instance; ///< mesh database
433  moab::Interface &moab = mb_instance; ///< mesh database interface
434  //! [Create MoAB]
435 
436  //! [Create MoFEM]
437  MoFEM::Core core(moab); ///< finite element database
438  MoFEM::Interface &m_field = core; ///< finite element database insterface
439  //! [Create MoFEM]
440 
441  //! [Example]
442  Example ex(m_field);
443  CHKERR ex.runProblem();
444  //! [Example]
445  }
446  CATCH_ERRORS;
447 
449 }
const double T
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:322
boost::shared_ptr< MatrixDouble > matInertiaPtr
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForce
structure for User Loop Methods on finite elementsIt can be used to calculate stiffness matrices,...
Volume finite element baseUser is implementing own operator at Gauss point level, by class derived fr...
Deprecated interface functions.
Problem manager is used to build and partition problems.
MoFEMErrorCode createCommonData()
[Set up problem]
MoFEMErrorCode setupProblem()
[Set up problem]
FTensor::Index< 'i', 3 > i
constexpr int SPACE_DIM
FTensor::Index< 'k', 3 > k
FTensor::Index< 'l', 3 > l
CoreTmp< 0 > Core
Definition: Core.hpp:1129
boost::shared_ptr< PostProcEle > postProc
ElementsAndOps< SPACE_DIM >::PostProcEle PostProcEle
Definition: contact.cpp:59
[Example]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:445
static double * ts_aa_ptr
#define MOFEM_LOG_TAG(channel, tag)
Tag channelTag channel tag is set until MOFEM_LOG_CHANNEL is called, then new tag can be set.
Definition: LogManager.hpp:334
constexpr double young_modulus
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
static constexpr int approx_order
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:303
MatrixDouble invJac
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:306
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:60
Approximate field valuse for given petsc vector.
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: helmholtz.cpp:113
Example(MoFEM::Interface &m_field)
[Push operators to pipeline]
MoFEMErrorCode outputResults()
[Solve]
Definition: helmholtz.cpp:256
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
Definition: AuxPETSc.hpp:238
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
Definition: helmholtz.cpp:236
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForce
Definition: elastic.cpp:39
static double * ts_time_ptr
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:278
MoFEMErrorCode postProcess()
function is run at the end of loop
DataForcesAndSourcesCore::EntData EntData
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:939
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:363
Kronecker Delta class symmetric.
MoFEMErrorCode runProblem()
[Run problem]
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
constexpr int order
MoFEMErrorCode checkResults()
[Postprocess results]
Calculate inverse of jacobian for face element.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
Post processing.
FTensor::Index< 'j', 3 > j
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: contact.cpp:95
PipelineManager interface.
#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.
#define CHKERR
Inline error check.
Definition: definitions.h:604
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
static char help[]
[Check]
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMMoFEM.cpp:915
constexpr double omega
Transform local reference derivatives of shape functions to global derivatives.
constexpr double rho
assemble the stiffness matrix
Postprocess on face.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpK
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
SmartPetscObj< DM > dM
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
Data on single entity (This is passed as argument to DataOperator::doWork)
MoFEMErrorCode readMesh()
[Read mesh]
Definition: helmholtz.cpp:80
MoFEMErrorCode assembleSystem()
[Applying essential BC]
Definition: helmholtz.cpp:161
Monitor(SmartPetscObj< DM > dm, boost::shared_ptr< PostProcEle > post_proc)
Core (interface) class.
Definition: Core.hpp:77
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:441
constexpr double poisson_ratio
continuous field
Definition: definitions.h:177
[Class definition]
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
std::string param_file
constexpr double bulk_modulus_K
constexpr double shear_modulus_G
boost::shared_ptr< MatrixDouble > matAccelerationPtr
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
Definition: contact.cpp:97
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
Face finite elementUser is implementing own operator at Gauss point level, by own object derived from...
int main(int argc, char *argv[])