v0.13.1
UnsaturatedFlow.hpp
Go to the documentation of this file.
1/** \file UnsaturatedFlow.hpp
2 * \brief Mix implementation of transport element
3 * \example UnsaturatedFlow.hpp
4 *
5 * \ingroup mofem_mix_transport_elem
6 *
7 */
8
9#ifndef __UNSATURATD_FLOW_HPP__
10#define __UNSATURATD_FLOW_HPP__
11
12namespace MixTransport {
13
14using PostProcEle = PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore>;
15
16/**
17 * \brief Generic material model for unsaturated water transport
18 *
19 * \note This is abstact class, no physical material model is implemented
20 * here.
21 */
23
24 virtual ~GenericMaterial() {}
25
26 static double ePsilon0; ///< Regularization parameter
27 static double ePsilon1; ///< Regularization parameter
28 static double scaleZ; ///< Scale z direction
29
30 double sCale; ///< Scale time dependent eq.
31
32 double h; ///< hydraulic head
33 double h_t; ///< rate of hydraulic head
34 // double h_flux_residual; ///< residual at point
35 double K; ///< Hydraulic conductivity [L/s]
36 double diffK; ///< Derivative of hydraulic conductivity [L/s * L^2/F]
37 double C; ///< Capacity [S^2/L^2]
38 double diffC; ///< Derivative of capacity [S^2/L^2 * L^2/F ]
39 double tHeta; ///< Water content
40 double Se; ///< Effective saturation
41
42 Range tEts; ///< Elements with this material
43
44 double x, y, z; ///< in meters (L)
45
46 /**
47 * \brief Initialize head
48 * @return value of head
49 */
50 virtual double initialPcEval() const = 0;
51 virtual void printMatParameters(const int id,
52 const std::string &prefix) const = 0;
53
54 virtual MoFEMErrorCode calK() {
56 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
57 "Not implemented how to calculate hydraulic conductivity");
59 }
60
63 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
64 "Not implemented how to calculate derivative of hydraulic "
65 "conductivity");
67 }
68
69 virtual MoFEMErrorCode calC() {
71 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
72 "Not implemented how to calculate capacity");
74 }
75
78 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
79 "Not implemented how to calculate capacity");
81 }
82
85 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
86 "Not implemented how to calculate capacity");
88 }
89
92 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
93 "Not implemented how to calculate capacity");
95 }
96};
97
98/**
99 * \brief Implementation of operators, problem and finite elements for
100 * unsaturated flow
101 */
103
104 DM dM; ///< Discrete manager for unsaturated flow problem
105
107 : MixTransportElement(m_field), dM(PETSC_NULL), lastEvalBcValEnt(0),
110
112 if (dM != PETSC_NULL) {
113 CHKERR DMDestroy(&dM);
114 CHKERRABORT(PETSC_COMM_WORLD, ierr);
115 }
116 }
117
118 typedef std::map<int, boost::shared_ptr<GenericMaterial>> MaterialsDoubleMap;
119 MaterialsDoubleMap dMatMap; ///< materials database
120
121 /**
122 * \brief For given element handle get material block Id
123 * @param ent finite element entity handle
124 * @param block_id reference to returned block id
125 * @return error code
126 */
127 virtual MoFEMErrorCode getMaterial(const EntityHandle ent,
128 int &block_id) const {
130 for (MaterialsDoubleMap::const_iterator mit = dMatMap.begin();
131 mit != dMatMap.end(); mit++) {
132 if (mit->second->tEts.find(ent) != mit->second->tEts.end()) {
133 block_id = mit->first;
135 }
136 }
138 "Element not found, no material data");
140 }
141
142 /**
143 * \brief Class storing information about boundary condition
144 */
145 struct BcData {
147 double fixValue;
148 boost::function<double(const double x, const double y, const double z)>
150 BcData() : hookFun(NULL) {}
151 };
152 typedef map<int, boost::shared_ptr<BcData>> BcMap;
153 BcMap bcValueMap; ///< Store boundary condition for head capillary pressure
154
155 EntityHandle lastEvalBcValEnt;
157
158 /**
159 * \brief Get value on boundary
160 * @param ent entity handle
161 * @param gg number of integration point
162 * @param x x-coordinate
163 * @param y y-coordinate
164 * @param z z-coordinate
165 * @param value returned value
166 * @return error code
167 */
168 MoFEMErrorCode getBcOnValues(const EntityHandle ent, const int gg,
169 const double x, const double y, const double z,
170 double &value) {
172 int block_id = -1;
173 if (lastEvalBcValEnt == ent) {
174 block_id = lastEvalBcBlockValId;
175 } else {
176 for (BcMap::iterator it = bcValueMap.begin(); it != bcValueMap.end();
177 it++) {
178 if (it->second->eNts.find(ent) != it->second->eNts.end()) {
179 block_id = it->first;
180 }
181 }
182 lastEvalBcValEnt = ent;
183 lastEvalBcBlockValId = block_id;
184 }
185 if (block_id >= 0) {
186 if (bcValueMap.at(block_id)->hookFun) {
187 value = bcValueMap.at(block_id)->hookFun(x, y, z);
188 } else {
189 value = bcValueMap.at(block_id)->fixValue;
190 }
191 } else {
192 value = 0;
193 }
195 }
196
198 EntityHandle lastEvalBcFluxEnt;
200
201 /**
202 * \brief essential (Neumann) boundary condition (set fluxes)
203 * @param ent handle to finite element entity
204 * @param x coord
205 * @param y coord
206 * @param z coord
207 * @param flux reference to flux which is set by function
208 * @return [description]
209 */
210 MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x,
211 const double y, const double z, double &flux) {
213 int block_id = -1;
214 if (lastEvalBcFluxEnt == ent) {
215 block_id = lastEvalBcBlockFluxId;
216 } else {
217 for (BcMap::iterator it = bcFluxMap.begin(); it != bcFluxMap.end();
218 it++) {
219 if (it->second->eNts.find(ent) != it->second->eNts.end()) {
220 block_id = it->first;
221 }
222 }
223 lastEvalBcFluxEnt = ent;
224 lastEvalBcBlockFluxId = block_id;
225 }
226 if (block_id >= 0) {
227 if (bcFluxMap.at(block_id)->hookFun) {
228 flux = bcFluxMap.at(block_id)->hookFun(x, y, z);
229 } else {
230 flux = bcFluxMap.at(block_id)->fixValue;
231 }
232 } else {
233 flux = 0;
234 }
236 }
237
238 /**
239 * \brief Evaluate boundary condition at the boundary
240 */
243
245 boost::shared_ptr<MethodForForceScaling> valueScale;
246
247 /**
248 * \brief Constructor
249 */
250 OpRhsBcOnValues(UnsaturatedFlowElement &ctx, const std::string fluxes_name,
251 boost::shared_ptr<MethodForForceScaling> &value_scale)
253 fluxes_name, UserDataOperator::OPROW),
254 cTx(ctx), valueScale(value_scale) {}
255
256 VectorDouble nF; ///< Vector of residuals
257
258 /**
259 * \brief Integrate boundary condition
260 * @param side local index of entity
261 * @param type type of entity
262 * @param data data on entity
263 * @return error code
264 */
268 if (data.getFieldData().size() == 0)
270 // Get EntityHandle of the finite element
271 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
272 // Resize and clear vector
273 nF.resize(data.getIndices().size());
274 nF.clear();
275 // Get number of integration points
276 int nb_gauss_pts = data.getN().size1();
277 for (int gg = 0; gg < nb_gauss_pts; gg++) {
278 double x, y, z;
279 x = getCoordsAtGaussPts()(gg, 0);
280 y = getCoordsAtGaussPts()(gg, 1);
281 z = getCoordsAtGaussPts()(gg, 2);
282 double value;
283 // get value of boundary condition
284 CHKERR cTx.getBcOnValues(fe_ent, gg, x, y, z, value);
285 const double w = getGaussPts()(2, gg) * 0.5;
286 const double beta = w * (value - z);
287 noalias(nF) += beta * prod(data.getVectorN<3>(gg), getNormal());
288 }
289 // Scale vector if history evaluating method is given
290 Vec f = getFEMethod()->ts_F;
291 if (valueScale) {
292 CHKERR valueScale->scaleNf(getFEMethod(), nF);
293 }
294 // Assemble vector
295 CHKERR VecSetValues(f, data.getIndices().size(), &data.getIndices()[0],
296 &nF[0], ADD_VALUES);
298 }
299 };
300
301 /**
302 * \brief Assemble flux residual
303 */
306
308
309 OpResidualFlux(UnsaturatedFlowElement &ctx, const std::string &flux_name)
310 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
311 flux_name, UserDataOperator::OPROW),
312 cTx(ctx) {}
313
316
320 const int nb_dofs = data.getIndices().size();
321 if (nb_dofs == 0)
323 nF.resize(nb_dofs, false);
324 nF.clear();
325 // Get EntityHandle of the finite element
326 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
327 // Get material block id
328 int block_id;
329 CHKERR cTx.getMaterial(fe_ent, block_id);
330 // Get material block
331 boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
332 // block_data->printMatParameters(block_id,"Read material");
333 // Get base function
334 auto t_n_hdiv = data.getFTensor1N<3>();
335 // Get pressure
337 // Get flux
338 auto t_flux = getFTensor1FromMat<3>(cTx.fluxesAtGaussPts);
339 // Coords at integration points
340 auto t_coords = getFTensor1CoordsAtGaussPts();
341 // Get integration weight
342 auto t_w = getFTensor0IntegrationWeight();
343 // Get volume
344 double vol = getVolume();
345
347 auto t_base_diff_hdiv = data.getFTensor2DiffN<3, 3>();
348 divVec.resize(nb_dofs, false);
349
350 // Get material parameters
351 int nb_gauss_pts = data.getN().size1();
352 for (int gg = 0; gg != nb_gauss_pts; gg++) {
353 // Get divergence
354 for (auto &v : divVec) {
355 v = t_base_diff_hdiv(i, i);
356 ++t_base_diff_hdiv;
357 }
358
359 const double alpha = t_w * vol;
360 block_data->h = t_h;
361 block_data->x = t_coords(0);
362 block_data->y = t_coords(1);
363 block_data->z = t_coords(2);
364 CHKERR block_data->calK();
365 const double K = block_data->K;
366 const double scaleZ = block_data->scaleZ;
367 const double z = t_coords(2); /// z-coordinate at Gauss pt
368 // Calculate pressure gradient
369 noalias(nF) -= alpha * (t_h - z * scaleZ) * divVec;
370 // Calculate presure gradient from flux
371 FTensor::Tensor0<double *> t_nf(&*nF.begin());
372 for (int rr = 0; rr != nb_dofs; rr++) {
373 t_nf += alpha * (1 / K) * (t_n_hdiv(i) * t_flux(i));
374 ++t_n_hdiv; // move to next base function
375 ++t_nf; // move to next element in vector
376 }
377 ++t_h; // move to next integration point
378 ++t_flux;
379 ++t_coords;
380 ++t_w;
381 }
382 // Assemble residual
383 CHKERR VecSetValues(getFEMethod()->ts_F, nb_dofs,
384 &*data.getIndices().begin(), &*nF.begin(),
385 ADD_VALUES);
387 }
388 };
389
390 /**
391 * Assemble mass residual
392 */
395
397
398 OpResidualMass(UnsaturatedFlowElement &ctx, const std::string &val_name)
399 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
400 val_name, UserDataOperator::OPROW),
401 cTx(ctx) {}
402
404
408 const int nb_dofs = data.getIndices().size();
409 if (nb_dofs == 0)
411 // Resize local element vector
412 nF.resize(nb_dofs, false);
413 nF.clear();
414 // Get EntityHandle of the finite element
415 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
416 // Get material block id
417 int block_id;
418 CHKERR cTx.getMaterial(fe_ent, block_id);
419 // Get material block
420 boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
421 // Get pressure
423 // Get pressure rate
425 // Flux divergence
427 // Get integration weight
428 auto t_w = getFTensor0IntegrationWeight();
429 // Coords at integration points
430 auto t_coords = getFTensor1CoordsAtGaussPts();
431 // Scale eq.
432 const double scale = block_data->sCale;
433 // Get volume
434 const double vol = getVolume();
435 // Get number of integration points
436 int nb_gauss_pts = data.getN().size1();
437 for (int gg = 0; gg != nb_gauss_pts; gg++) {
438 const double alpha = t_w * vol * scale;
439 block_data->h = t_h;
440 block_data->x = t_coords(0);
441 block_data->y = t_coords(1);
442 block_data->z = t_coords(2);
443 CHKERR block_data->calC();
444 const double C = block_data->C;
445 // Calculate flux conservation
446 noalias(nF) += (alpha * (t_div_flux + C * t_h_t)) * data.getN(gg);
447 ++t_h;
448 ++t_h_t;
449 ++t_div_flux;
450 ++t_coords;
451 ++t_w;
452 }
453 // Assemble local vector
454 Vec f = getFEMethod()->ts_F;
455 CHKERR VecSetValues(f, nb_dofs, &*data.getIndices().begin(), &*nF.begin(),
456 ADD_VALUES);
458 }
459 };
460
463
465
467 const std::string flux_name)
468 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
469 flux_name, flux_name, UserDataOperator::OPROWCOL),
470 cTx(ctx) {
471 sYmm = true;
472 }
473
475
477
478 /**
479 * \brief Assemble matrix
480 * @param row_side local index of row entity on element
481 * @param col_side local index of col entity on element
482 * @param row_type type of row entity, f.e. MBVERTEX, MBEDGE, or MBTET
483 * @param col_type type of col entity, f.e. MBVERTEX, MBEDGE, or MBTET
484 * @param row_data data for row
485 * @param col_data data for col
486 * @return error code
487 */
488 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
489 EntityType col_type,
491 EntitiesFieldData::EntData &col_data) {
493 const int nb_row = row_data.getIndices().size();
494 const int nb_col = col_data.getIndices().size();
495 if (nb_row == 0)
497 if (nb_col == 0)
499 nN.resize(nb_row, nb_col, false);
500 nN.clear();
501 // Get EntityHandle of the finite element
502 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
503 // Get material block id
504 int block_id;
505 CHKERR cTx.getMaterial(fe_ent, block_id);
506 // Get material block
507 boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
508 // Get pressure
510 // Coords at integration points
511 auto t_coords = getFTensor1CoordsAtGaussPts();
512 // Get base functions
513 auto t_n_hdiv_row = row_data.getFTensor1N<3>();
514 // Get integration weight
515 auto t_w = getFTensor0IntegrationWeight();
516 // Get volume
517 const double vol = getVolume();
518 int nb_gauss_pts = row_data.getN().size1();
519 for (int gg = 0; gg != nb_gauss_pts; gg++) {
520 block_data->h = t_h;
521 block_data->x = t_coords(0);
522 block_data->y = t_coords(1);
523 block_data->z = t_coords(2);
524 CHKERR block_data->calK();
525 const double K = block_data->K;
526 // get integration weight and multiply by element volume
527 const double alpha = t_w * vol;
528 const double beta = alpha * (1 / K);
529 FTensor::Tensor0<double *> t_a(&*nN.data().begin());
530 for (int kk = 0; kk != nb_row; kk++) {
531 auto t_n_hdiv_col = col_data.getFTensor1N<3>(gg, 0);
532 for (int ll = 0; ll != nb_col; ll++) {
533 t_a += beta * (t_n_hdiv_row(j) * t_n_hdiv_col(j));
534 ++t_n_hdiv_col;
535 ++t_a;
536 }
537 ++t_n_hdiv_row;
538 }
539 ++t_coords;
540 ++t_h;
541 ++t_w;
542 }
543 Mat a = getFEMethod()->ts_B;
544 CHKERR MatSetValues(a, nb_row, &*row_data.getIndices().begin(), nb_col,
545 &*col_data.getIndices().begin(), &*nN.data().begin(),
546 ADD_VALUES);
547 // matrix is symmetric, assemble other part
548 if (row_side != col_side || row_type != col_type) {
549 transNN.resize(nb_col, nb_row);
550 noalias(transNN) = trans(nN);
551 CHKERR MatSetValues(a, nb_col, &*col_data.getIndices().begin(), nb_row,
552 &*row_data.getIndices().begin(),
553 &*transNN.data().begin(), ADD_VALUES);
554 }
556 }
557 };
558
561
563
564 OpVU_L2L2(UnsaturatedFlowElement &ctx, const std::string value_name)
565 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
566 value_name, value_name, UserDataOperator::OPROWCOL),
567 cTx(ctx) {
568 sYmm = true;
569 }
570
572
573 /**
574 * \brief Assemble matrix
575 * @param row_side local index of row entity on element
576 * @param col_side local index of col entity on element
577 * @param row_type type of row entity, f.e. MBVERTEX, MBEDGE, or MBTET
578 * @param col_type type of col entity, f.e. MBVERTEX, MBEDGE, or MBTET
579 * @param row_data data for row
580 * @param col_data data for col
581 * @return error code
582 */
583 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
584 EntityType col_type,
586 EntitiesFieldData::EntData &col_data) {
588 int nb_row = row_data.getIndices().size();
589 int nb_col = col_data.getIndices().size();
590 if (nb_row == 0)
592 if (nb_col == 0)
594 nN.resize(nb_row, nb_col, false);
595 nN.clear();
596 // Get EntityHandle of the finite element
597 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
598 // Get material block id
599 int block_id;
600 CHKERR cTx.getMaterial(fe_ent, block_id);
601 // Get material block
602 boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
603 // Get pressure
605 // // Get pressure
606 // auto t_flux_residual = getFTensor0FromVec(*cTx.resAtGaussPts);
607 // Get pressure rate
609 // Get integration weight
610 auto t_w = getFTensor0IntegrationWeight();
611 // Coords at integration points
612 auto t_coords = getFTensor1CoordsAtGaussPts();
613 // Scale eq.
614 const double scale = block_data->sCale;
615 // Time step factor
616 double ts_a = getFEMethod()->ts_a;
617 // get volume
618 const double vol = getVolume();
619 int nb_gauss_pts = row_data.getN().size1();
620 // get base functions on rows
621 auto t_n_row = row_data.getFTensor0N();
622 for (int gg = 0; gg != nb_gauss_pts; gg++) {
623 // get integration weight and multiply by element volume
624 double alpha = t_w * vol * scale;
625 // evaluate material model at integration points
626 // to calculate capacity and tangent of capacity term
627 block_data->h = t_h;
628 block_data->h_t = t_h_t;
629 block_data->x = t_coords(0);
630 block_data->y = t_coords(1);
631 block_data->z = t_coords(2);
632 CHKERR block_data->calC();
633 CHKERR block_data->calDiffC();
634 const double C = block_data->C;
635 const double diffC = block_data->diffC;
636 // assemble local entity tangent matrix
638 &*nN.data().begin());
639 // iterate base functions on rows
640 for (int kk = 0; kk != nb_row; kk++) {
641 // get first base function on column at integration point gg
642 auto t_n_col = col_data.getFTensor0N(gg, 0);
643 // iterate base functions on columns
644 for (int ll = 0; ll != nb_col; ll++) {
645 // assemble elements of local matrix
646 t_a += (alpha * (C * ts_a + diffC * t_h_t)) * t_n_row * t_n_col;
647 ++t_n_col; // move to next base function on column
648 ++t_a; // move to next element in local tangent matrix
649 }
650 ++t_n_row; // move to next base function on row
651 }
652 ++t_w; // move to next integration weight
653 ++t_coords; // move to next coordinate at integration point
654 ++t_h; // move to next capillary head at integration point
655 // ++t_flux_residual;
656 ++t_h_t; // move to next capillary head rate at integration point
657 }
658 Mat a = getFEMethod()->ts_B;
659 CHKERR MatSetValues(a, nb_row, &row_data.getIndices()[0], nb_col,
660 &col_data.getIndices()[0], &*nN.data().begin(),
661 ADD_VALUES);
663 }
664 };
665
668
670
671 /**
672 * \brief Constructor
673 */
675 const std::string &val_name_row,
676 const std::string &flux_name_col)
677 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
678 val_name_row, flux_name_col, UserDataOperator::OPROWCOL, false),
679 cTx(ctx) {}
680
683
684 /**
685 * \brief Do calculations
686 * @param row_side local index of entity on row
687 * @param col_side local index of entity on column
688 * @param row_type type of row entity
689 * @param col_type type of col entity
690 * @param row_data row data structure carrying information about base
691 * functions, DOFs indices, etc.
692 * @param col_data column data structure carrying information about base
693 * functions, DOFs indices, etc.
694 * @return error code
695 */
696 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
697 EntityType col_type,
699 EntitiesFieldData::EntData &col_data) {
701 int nb_row = row_data.getFieldData().size();
702 int nb_col = col_data.getFieldData().size();
703 if (nb_row == 0)
705 if (nb_col == 0)
707 // Get EntityHandle of the finite element
708 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
709 // Get material block id
710 int block_id;
711 CHKERR cTx.getMaterial(fe_ent, block_id);
712 // Get material block
713 boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
714 nN.resize(nb_row, nb_col, false);
715 divVec.resize(nb_col, false);
716 nN.clear();
717 // take derivatives of base functions
719 auto t_base_diff_hdiv = col_data.getFTensor2DiffN<3, 3>();
720 // Scale eq.
721 const double scale = block_data->sCale;
722 int nb_gauss_pts = row_data.getN().size1();
723 for (int gg = 0; gg < nb_gauss_pts; gg++) {
724 double alpha = getGaussPts()(3, gg) * getVolume() * scale;
725 for (auto &v : divVec) {
726 v = t_base_diff_hdiv(i, i);
727 ++t_base_diff_hdiv;
728 }
729 noalias(nN) += alpha * outer_prod(row_data.getN(gg), divVec);
730 }
731 CHKERR MatSetValues(getFEMethod()->ts_B, nb_row,
732 &row_data.getIndices()[0], nb_col,
733 &col_data.getIndices()[0], &nN(0, 0), ADD_VALUES);
735 }
736 };
737
740
742
743 /**
744 * \brief Constructor
745 */
747 const std::string &flux_name_col,
748 const std::string &val_name_row)
749 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
750 flux_name_col, val_name_row, UserDataOperator::OPROWCOL, false),
751 cTx(ctx) {}
752
756
757 /**
758 * \brief Do calculations
759 * @param row_side local index of entity on row
760 * @param col_side local index of entity on column
761 * @param row_type type of row entity
762 * @param col_type type of col entity
763 * @param row_data row data structure carrying information about base
764 * functions, DOFs indices, etc.
765 * @param col_data column data structure carrying information about base
766 * functions, DOFs indices, etc.
767 * @return error code
768 */
769 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
770 EntityType col_type,
772 EntitiesFieldData::EntData &col_data) {
774 int nb_row = row_data.getFieldData().size();
775 int nb_col = col_data.getFieldData().size();
776 if (nb_row == 0)
778 if (nb_col == 0)
780 nN.resize(nb_row, nb_col, false);
781 divVec.resize(nb_row, false);
782 nN.clear();
783 // Get EntityHandle of the finite element
784 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
785 // Get material block id
786 int block_id;
787 CHKERR cTx.getMaterial(fe_ent, block_id);
788 // Get material block
789 boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
790 // Get pressure
792 // Get flux
793 auto t_flux = getFTensor1FromMat<3>(cTx.fluxesAtGaussPts);
794 // Coords at integration points
795 auto t_coords = getFTensor1CoordsAtGaussPts();
796 // Get integration weight
797 auto t_w = getFTensor0IntegrationWeight();
798 // Get base function
799 auto t_n_hdiv_row = row_data.getFTensor1N<3>();
800 // Get derivatives of base functions
802 auto t_base_diff_hdiv = row_data.getFTensor2DiffN<3, 3>();
803 // Get volume
804 double vol = getVolume();
805 int nb_gauss_pts = row_data.getN().size1();
806 for (int gg = 0; gg != nb_gauss_pts; gg++) {
807 block_data->h = t_h;
808 block_data->x = t_coords(0);
809 block_data->y = t_coords(1);
810 block_data->z = t_coords(2);
811 CHKERR block_data->calK();
812 CHKERR block_data->calDiffK();
813 const double K = block_data->K;
814 // const double z = t_coords(2);
815 const double KK = K * K;
816 const double diffK = block_data->diffK;
817 double alpha = t_w * vol;
818 for (auto &v : divVec) {
819 v = t_base_diff_hdiv(i, i);
820 ++t_base_diff_hdiv;
821 }
822 noalias(nN) -= alpha * outer_prod(divVec, col_data.getN(gg));
823 FTensor::Tensor0<double *> t_a(&*nN.data().begin());
824 for (int rr = 0; rr != nb_row; rr++) {
825 double beta = alpha * (-diffK / KK) * (t_n_hdiv_row(i) * t_flux(i));
826 auto t_n_col = col_data.getFTensor0N(gg, 0);
827 for (int cc = 0; cc != nb_col; cc++) {
828 t_a += beta * t_n_col;
829 ++t_n_col;
830 ++t_a;
831 }
832 ++t_n_hdiv_row;
833 }
834 ++t_w;
835 ++t_coords;
836 ++t_h;
837 ++t_flux;
838 }
839 CHKERR MatSetValues(getFEMethod()->ts_B, nb_row,
840 &row_data.getIndices()[0], nb_col,
841 &col_data.getIndices()[0], &nN(0, 0), ADD_VALUES);
843 }
844 };
845
850 const std::string &val_name)
851 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
852 val_name, UserDataOperator::OPROW),
853 cTx(ctx) {}
854
857
861 if (data.getFieldData().size() == 0)
863 int nb_dofs = data.getFieldData().size();
864 int nb_gauss_pts = data.getN().size1();
865 if (nb_dofs != static_cast<int>(data.getN().size2())) {
866 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
867 "wrong number of dofs");
868 }
869 nN.resize(nb_dofs, nb_dofs, false);
870 nF.resize(nb_dofs, false);
871 nN.clear();
872 nF.clear();
873
874 // Get EntityHandle of the finite element
875 EntityHandle fe_ent = getFEEntityHandle();
876 // Get material block id
877 int block_id;
878 CHKERR cTx.getMaterial(fe_ent, block_id);
879 // Get material block
880 boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
881
882 // loop over integration points
883 for (int gg = 0; gg < nb_gauss_pts; gg++) {
884 // get coordinates at integration point
885 block_data->x = getCoordsAtGaussPts()(gg, 0);
886 block_data->y = getCoordsAtGaussPts()(gg, 1);
887 block_data->z = getCoordsAtGaussPts()(gg, 2);
888 // get weight for integration rule
889 double alpha = getGaussPts()(2, gg) * getVolume();
890 nN += alpha * outer_prod(data.getN(gg), data.getN(gg));
891 nF += alpha * block_data->initialPcEval() * data.getN(gg);
892 }
893
894 // factor matrix
896 // solve local problem
897 cholesky_solve(nN, nF, ublas::lower());
898
899 // set solution to vector
900 CHKERR VecSetValues(cTx.D1, nb_dofs, &*data.getIndices().begin(),
901 &*nF.begin(), INSERT_VALUES);
902
904 }
905 };
906
910
911 /**
912 * \brief Constructor
913 */
915 const std::string fluxes_name)
917 fluxes_name, UserDataOperator::OPROW),
918 cTx(ctx) {}
919
921
922 /**
923 * \brief Integrate boundary condition
924 * @param side local index of entity
925 * @param type type of entity
926 * @param data data on entity
927 * @return error code
928 */
932 int nb_dofs = data.getFieldData().size();
933 if (nb_dofs == 0)
935 // Get base function
936 auto t_n_hdiv = data.getFTensor1N<3>();
937 // get normal of face
938 auto t_normal = getFTensor1Normal();
939 // Integration weight
940 auto t_w = getFTensor0IntegrationWeight();
941 double flux_on_entity = 0;
942 int nb_gauss_pts = data.getN().size1();
943 for (int gg = 0; gg < nb_gauss_pts; gg++) {
944 auto t_data = data.getFTensor0FieldData();
945 for (int rr = 0; rr != nb_dofs; rr++) {
946 flux_on_entity -= (0.5 * t_data * t_w) * (t_n_hdiv(i) * t_normal(i));
947 ++t_n_hdiv;
948 ++t_data;
949 }
950 ++t_w;
951 }
952 CHKERR VecSetValue(cTx.ghostFlux, 0, flux_on_entity, ADD_VALUES);
954 }
955 };
956
957 /**
958 * Operator used to post-process results for unsaturated infiltration problem.
959 * Operator should with element for post-processing results, i.e.
960 * PostProcVolumeOnRefinedMesh
961 */
964
967 std::vector<EntityHandle> &mapGaussPts;
968
970 moab::Interface &post_proc_mesh,
971 std::vector<EntityHandle> &map_gauss_pts,
972 const std::string field_name)
973 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
975 cTx(ctx), postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts) {}
976
980 int nb_dofs = data.getFieldData().size();
981 if (nb_dofs == 0)
983
984 // Get EntityHandle of the finite element
985 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
986 // Get material block id
987 int block_id;
988 CHKERR cTx.getMaterial(fe_ent, block_id);
989 // Get material block
990 boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
991
992 // Set bloc Id
993 Tag th_id;
994 int def_block_id = -1;
995 CHKERR postProcMesh.tag_get_handle("BLOCK_ID", 1, MB_TYPE_INTEGER, th_id,
996 MB_TAG_CREAT | MB_TAG_SPARSE,
997 &def_block_id);
998
999 // Create mesh tag. Tags are created on post-processing mesh and
1000 // visable in post-processor, e.g. Paraview
1001 double zero = 0;
1002 Tag th_theta;
1003 CHKERR postProcMesh.tag_get_handle("THETA", 1, MB_TYPE_DOUBLE, th_theta,
1004 MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
1005 Tag th_se;
1006 CHKERR postProcMesh.tag_get_handle("Se", 1, MB_TYPE_DOUBLE, th_se,
1007 MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
1008 // Tag th_ks;
1009 // CHKERR postProcMesh.tag_get_handle(
1010 // "Ks",1,MB_TYPE_DOUBLE,th_ks,
1011 // MB_TAG_CREAT|MB_TAG_SPARSE,&zero
1012 // );
1013 // CHKERR postProcMesh.tag_set_data(th_ks,&fe_ent,1,&block_data->Ks);
1014 Tag th_k;
1015 CHKERR postProcMesh.tag_get_handle("K", 1, MB_TYPE_DOUBLE, th_k,
1016 MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
1017 Tag th_c;
1018 CHKERR postProcMesh.tag_get_handle("C", 1, MB_TYPE_DOUBLE, th_c,
1019 MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
1020
1021 // Get pressure at integration points
1023 // Coords at integration points
1024 auto t_coords = getFTensor1CoordsAtGaussPts();
1025
1026 int nb_gauss_pts = data.getN().size1();
1027 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1028 block_data->h = t_h;
1029 block_data->x = t_coords(0);
1030 block_data->y = t_coords(1);
1031 block_data->z = t_coords(2);
1032 // Calculate theta (water content) and save it on mesh tags
1033 CHKERR block_data->calTheta();
1034 double theta = block_data->tHeta;
1035 CHKERR postProcMesh.tag_set_data(th_theta, &mapGaussPts[gg], 1, &theta);
1036 CHKERR block_data->calSe();
1037 // Calculate Se (effective saturation and save it on the mesh tags)
1038 double Se = block_data->Se;
1039 CHKERR postProcMesh.tag_set_data(th_se, &mapGaussPts[gg], 1, &Se);
1040 // Calculate K (hydraulic conductivity) and save it on the mesh tags
1041 CHKERR block_data->calK();
1042 double K = block_data->K;
1043 CHKERR postProcMesh.tag_set_data(th_k, &mapGaussPts[gg], 1, &K);
1044 // Calculate water capacity and save it on the mesh tags
1045 CHKERR block_data->calC();
1046 double C = block_data->C;
1047 CHKERR postProcMesh.tag_set_data(th_c, &mapGaussPts[gg], 1, &C);
1048
1049 // Set block iD
1050 CHKERR postProcMesh.tag_set_data(th_id, &mapGaussPts[gg], 1, &block_id);
1051
1052 ++t_h;
1053 ++t_coords;
1054 }
1055
1057 }
1058 };
1059
1060 /**
1061 * Finite element implementation called by TS monitor. Element calls other
1062 * finite elements to evaluate material properties and save results on the
1063 * mesh.
1064 *
1065 * \note Element overloaded only FEMethod::postProcess methods where other
1066 * elements are called.
1067 */
1069
1071 boost::shared_ptr<PostProcEle> postProc;
1072 boost::shared_ptr<ForcesAndSourcesCore> fluxIntegrate;
1073
1074 const int fRequency;
1075
1077 boost::shared_ptr<PostProcEle> &post_proc,
1078 boost::shared_ptr<ForcesAndSourcesCore> flux_Integrate,
1079 const int frequency)
1080 : cTx(ctx), postProc(post_proc), fluxIntegrate(flux_Integrate),
1081 fRequency(frequency) {}
1082
1086 }
1087
1091 }
1092
1095
1096 // Get time step
1097 int step;
1098 CHKERR TSGetTimeStepNumber(ts, &step);
1099
1100 if ((step) % fRequency == 0) {
1101 // Post-process results and save in the file
1102 PetscPrintf(PETSC_COMM_WORLD, "Output results %d - %d\n", step,
1103 fRequency);
1105 CHKERR postProc->writeFile(
1106 string("out_") + boost::lexical_cast<std::string>(step) + ".h5m");
1107 }
1108
1109 // Integrate fluxes on faces where pressure head is applied
1110 CHKERR VecZeroEntries(cTx.ghostFlux);
1111 CHKERR VecGhostUpdateBegin(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
1112 CHKERR VecGhostUpdateEnd(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
1113 // Run finite element to integrate fluxes
1115 CHKERR VecAssemblyBegin(cTx.ghostFlux);
1116 CHKERR VecAssemblyEnd(cTx.ghostFlux);
1117 // accumulate errors from processors
1118 CHKERR VecGhostUpdateBegin(cTx.ghostFlux, ADD_VALUES, SCATTER_REVERSE);
1119 CHKERR VecGhostUpdateEnd(cTx.ghostFlux, ADD_VALUES, SCATTER_REVERSE);
1120 // scatter errors to all processors
1121 CHKERR VecGhostUpdateBegin(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
1122 CHKERR VecGhostUpdateEnd(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
1123 double *ghost_flux;
1124 CHKERR VecGetArray(cTx.ghostFlux, &ghost_flux);
1125 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Flux at time %6.4g %6.4g\n", ts_t,
1126 ghost_flux[0]);
1127 CHKERR VecRestoreArray(cTx.ghostFlux, &ghost_flux);
1128
1130 }
1131 };
1132
1133 /// \brief add fields
1134 MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes,
1135 const int order) {
1137 // Fields
1140 CHKERR mField.add_field(values + "_t", L2, AINSWORTH_LEGENDRE_BASE, 1);
1141 // CHKERR mField.add_field(fluxes+"_residual",L2,AINSWORTH_LEGENDRE_BASE,1);
1142
1143 // meshset consisting all entities in mesh
1144 EntityHandle root_set = mField.get_moab().get_root_set();
1145 // add entities to field
1146
1148 if (it->getName().compare(0, 4, "SOIL") != 0)
1149 continue;
1150 CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
1151 MBTET, fluxes);
1152 CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
1153 MBTET, values);
1154 CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
1155 MBTET, values + "_t");
1156 // CHKERR mField.add_ents_to_field_by_type(
1157 // dMatMap[it->getMeshsetId()]->tEts,MBTET,fluxes+"_residual"
1158 // );
1159 }
1160
1161 CHKERR mField.set_field_order(root_set, MBTET, fluxes, order + 1);
1162 CHKERR mField.set_field_order(root_set, MBTRI, fluxes, order + 1);
1163 CHKERR mField.set_field_order(root_set, MBTET, values, order);
1164 CHKERR mField.set_field_order(root_set, MBTET, values + "_t", order);
1165 // CHKERR mField.set_field_order(root_set,MBTET,fluxes+"_residual",order);
1167 }
1168
1169 /// \brief add finite elements
1170 MoFEMErrorCode addFiniteElements(const std::string &fluxes_name,
1171 const std::string &values_name) {
1173
1174 // Define element "MIX". Note that this element will work with fluxes_name
1175 // and values_name. This reflect bilinear form for the problem
1184 values_name + "_t");
1185 // CHKERR
1186 // mField.modify_finite_element_add_field_data("MIX",fluxes_name+"_residual");
1187
1189 if (it->getName().compare(0, 4, "SOIL") != 0)
1190 continue;
1192 dMatMap[it->getMeshsetId()]->tEts, MBTET, "MIX");
1193 }
1194
1195 // Define element to integrate natural boundary conditions, i.e. set values.
1196 CHKERR mField.add_finite_element("MIX_BCVALUE", MF_ZERO);
1198 fluxes_name);
1200 fluxes_name);
1202 fluxes_name);
1204 values_name);
1205
1207 if (it->getName().compare(0, 4, "HEAD") != 0)
1208 continue;
1210 bcValueMap[it->getMeshsetId()]->eNts, MBTRI, "MIX_BCVALUE");
1211 }
1212
1213 // Define element to apply essential boundary conditions.
1214 CHKERR mField.add_finite_element("MIX_BCFLUX", MF_ZERO);
1216 fluxes_name);
1218 fluxes_name);
1220 fluxes_name);
1222 values_name);
1223
1225 if (it->getName().compare(0, 4, "FLUX") != 0)
1226 continue;
1228 bcFluxMap[it->getMeshsetId()]->eNts, MBTRI, "MIX_BCFLUX");
1229 }
1230
1232 }
1233
1234 /**
1235 * \brief Build problem
1236 * @param ref_level mesh refinement on which mesh problem you like to built.
1237 * @return error code
1238 */
1240 BitRefLevel ref_level = BitRefLevel().set(0)) {
1242
1243 // Build fields
1245 // Build finite elements
1247 CHKERR mField.build_finite_elements("MIX_BCFLUX");
1248 CHKERR mField.build_finite_elements("MIX_BCVALUE");
1249 // Build adjacencies of degrees of freedom and elements
1250 CHKERR mField.build_adjacencies(ref_level);
1251
1252 // create DM instance
1253 CHKERR DMCreate(PETSC_COMM_WORLD, &dM);
1254 // setting that DM is type of DMMOFEM, i.e. MOFEM implementation manages DM
1255 CHKERR DMSetType(dM, "DMMOFEM");
1256 // mesh is portioned, each process keeps only part of problem
1257 CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
1258 // creates problem in DM
1259 CHKERR DMMoFEMCreateMoFEM(dM, &mField, "MIX", ref_level);
1260 // discretised problem creates square matrix (that makes some optimizations)
1261 CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
1262 // set DM options from command line
1263 CHKERR DMSetFromOptions(dM);
1264 // add finite elements
1265 CHKERR DMMoFEMAddElement(dM, "MIX");
1266 CHKERR DMMoFEMAddElement(dM, "MIX_BCFLUX");
1267 CHKERR DMMoFEMAddElement(dM, "MIX_BCVALUE");
1268 // constructor data structures
1269 CHKERR DMSetUp(dM);
1270
1271 // remove zero flux dofs
1272 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
1273 "MIX", "FLUXES", zero_flux_ents);
1274
1275 PetscSection section;
1276 CHKERR mField.getInterface<ISManager>()->sectionCreate("MIX", &section);
1277 CHKERR DMSetSection(dM, section);
1278 // CHKERR PetscSectionView(section,PETSC_VIEWER_STDOUT_WORLD);
1279 CHKERR PetscSectionDestroy(&section);
1280
1282 }
1283
1284 boost::shared_ptr<ForcesAndSourcesCore>
1285 feFaceBc; ///< Elemnet to calculate essential bc
1286 boost::shared_ptr<ForcesAndSourcesCore>
1287 feFaceRhs; ///< Face element apply natural bc
1288 boost::shared_ptr<ForcesAndSourcesCore>
1289 feVolInitialPc; ///< Calculate inital boundary conditions
1290 boost::shared_ptr<ForcesAndSourcesCore>
1291 feVolRhs; ///< Assemble residual vector
1292 boost::shared_ptr<ForcesAndSourcesCore> feVolLhs; ///< Assemble tangent matrix
1293 boost::shared_ptr<MethodForForceScaling>
1294 scaleMethodFlux; ///< Method scaling fluxes
1295 boost::shared_ptr<MethodForForceScaling>
1296 scaleMethodValue; ///< Method scaling values
1297 boost::shared_ptr<FEMethod> tsMonitor; ///< Element used by TS monitor to
1298 ///< postprocess results at time step
1299
1300 boost::shared_ptr<VectorDouble>
1301 headRateAtGaussPts; ///< Vector keeps head rate
1302 // boost::shared_ptr<VectorDouble> resAtGaussPts; ///< Residual field
1303
1304 /**
1305 * \brief Set integration rule to volume elements
1306 *
1307 */
1308 struct VolRule {
1309 int operator()(int, int, int p_data) const { return 2 * p_data + p_data; }
1310 };
1311 /**
1312 * \brief Set integration rule to boundary elements
1313 *
1314 */
1315 struct FaceRule {
1316 int operator()(int p_row, int p_col, int p_data) const {
1317 return 2 * p_data;
1318 }
1319 };
1320
1321 std::vector<int> bcVecIds;
1323
1324 /**
1325 * \brief Pre-peprocessing
1326 * Set head pressute rate and get inital essential boundary conditions
1327 */
1330 boost::shared_ptr<ForcesAndSourcesCore> fePtr;
1331 // std::string mArk;
1332
1335 boost::shared_ptr<ForcesAndSourcesCore> &fe_ptr /*,std::string mark*/
1336 )
1337 : cTx(ctx), fePtr(fe_ptr) /*,mArk(mark)*/ {}
1340 // Update pressure rates
1341 CHKERR fePtr->mField.getInterface<VecManager>()->setOtherLocalGhostVector(
1342 fePtr->problemPtr, "VALUES", string("VALUES") + "_t", ROW,
1343 fePtr->ts_u_t, INSERT_VALUES, SCATTER_REVERSE);
1344 switch (fePtr->ts_ctx) {
1345 case TSMethod::CTX_TSSETIFUNCTION:
1346 CHKERR VecSetOption(fePtr->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
1347 PETSC_TRUE);
1348 if (!cTx.bcIndices.empty()) {
1349 double scale;
1350 CHKERR cTx.scaleMethodFlux->getForceScale(fePtr->ts_t, scale);
1351 if (cTx.bcVecIds.size() != cTx.bcIndices.size()) {
1352 cTx.bcVecIds.insert(cTx.bcVecIds.begin(), cTx.bcIndices.begin(),
1353 cTx.bcIndices.end());
1354 cTx.bcVecVals.resize(cTx.bcVecIds.size(), false);
1355 cTx.vecValsOnBc.resize(cTx.bcVecIds.size(), false);
1356 }
1357 CHKERR VecGetValues(cTx.D0, cTx.bcVecIds.size(),
1358 &*cTx.bcVecIds.begin(), &*cTx.bcVecVals.begin());
1359 CHKERR VecGetValues(fePtr->ts_u, cTx.bcVecIds.size(),
1360 &*cTx.bcVecIds.begin(),
1361 &*cTx.vecValsOnBc.begin());
1362 cTx.bcVecVals *= scale;
1363 VectorDouble::iterator vit = cTx.bcVecVals.begin();
1364 const NumeredDofEntity *dof_ptr;
1365 for (std::vector<int>::iterator it = cTx.bcVecIds.begin();
1366 it != cTx.bcVecIds.end(); it++, vit++) {
1367 if (auto dof_ptr =
1368 fePtr->problemPtr->getColDofsByPetscGlobalDofIdx(*it)
1369 .lock())
1370 dof_ptr->getFieldData() = *vit;
1371 }
1372 } else {
1373 cTx.bcVecIds.resize(0);
1374 cTx.bcVecVals.resize(0);
1375 cTx.vecValsOnBc.resize(0);
1376 }
1377 break;
1378 default:
1379 // don nothing
1380 break;
1381 }
1383 }
1384 };
1385
1386 /**
1387 * \brief Post proces method for volume element
1388 * Assemble vectors and matrices and apply essential boundary conditions
1389 */
1392 boost::shared_ptr<ForcesAndSourcesCore> fePtr;
1393 // std::string mArk;
1396 boost::shared_ptr<ForcesAndSourcesCore> &fe_ptr //,std::string mark
1397 )
1398 : cTx(ctx), fePtr(fe_ptr) /*,mArk(mark)*/ {}
1401 switch (fePtr->ts_ctx) {
1402 case TSMethod::CTX_TSSETIJACOBIAN: {
1403 CHKERR MatAssemblyBegin(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
1404 CHKERR MatAssemblyEnd(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
1405 // MatView(fePtr->ts_B,PETSC_VIEWER_DRAW_WORLD);
1406 // std::string wait;
1407 // std::cin >> wait;
1408 CHKERR MatZeroRowsColumns(fePtr->ts_B, cTx.bcVecIds.size(),
1409 &*cTx.bcVecIds.begin(), 1, PETSC_NULL,
1410 PETSC_NULL);
1411 CHKERR MatAssemblyBegin(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
1412 CHKERR MatAssemblyEnd(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
1413 // MatView(fePtr->ts_B,PETSC_VIEWER_DRAW_WORLD);
1414 // std::string wait;
1415 // std::cin >> wait;
1416 } break;
1417 case TSMethod::CTX_TSSETIFUNCTION: {
1418 CHKERR VecAssemblyBegin(fePtr->ts_F);
1419 CHKERR VecAssemblyEnd(fePtr->ts_F);
1420 if (!cTx.bcVecIds.empty()) {
1422 CHKERR VecSetValues(fePtr->ts_F, cTx.bcVecIds.size(),
1423 &*cTx.bcVecIds.begin(), &*cTx.vecValsOnBc.begin(),
1424 INSERT_VALUES);
1425 }
1426 CHKERR VecAssemblyBegin(fePtr->ts_F);
1427 CHKERR VecAssemblyEnd(fePtr->ts_F);
1428 } break;
1429 default:
1430 // don nothing
1431 break;
1432 }
1434 }
1435 };
1436
1437 /**
1438 * \brief Create finite element instances
1439 * @param vol_rule integration rule for volume element
1440 * @param face_rule integration rule for boundary element
1441 * @return error code
1442 */
1444 setFiniteElements(ForcesAndSourcesCore::RuleHookFun vol_rule = VolRule(),
1445 ForcesAndSourcesCore::RuleHookFun face_rule = FaceRule()) {
1447
1448 // create finite element instances
1449 feFaceBc = boost::shared_ptr<ForcesAndSourcesCore>(
1451 feFaceRhs = boost::shared_ptr<ForcesAndSourcesCore>(
1453 feVolInitialPc = boost::shared_ptr<ForcesAndSourcesCore>(
1454 new VolumeElementForcesAndSourcesCore(mField));
1455 feVolLhs = boost::shared_ptr<ForcesAndSourcesCore>(
1456 new VolumeElementForcesAndSourcesCore(mField));
1457 feVolRhs = boost::shared_ptr<ForcesAndSourcesCore>(
1458 new VolumeElementForcesAndSourcesCore(mField));
1459 // set integration rule to elements
1460 feFaceBc->getRuleHook = face_rule;
1461 feFaceRhs->getRuleHook = face_rule;
1462 feVolInitialPc->getRuleHook = vol_rule;
1463 feVolLhs->getRuleHook = vol_rule;
1464 feVolRhs->getRuleHook = vol_rule;
1465 // set function hook for finite element postprocessing stage
1466 feVolRhs->preProcessHook = preProcessVol(*this, feVolRhs);
1467 feVolLhs->preProcessHook = preProcessVol(*this, feVolLhs);
1468 feVolRhs->postProcessHook = postProcessVol(*this, feVolRhs);
1469 feVolLhs->postProcessHook = postProcessVol(*this, feVolLhs);
1470
1471 // Set Piola transform
1472 feFaceBc->getOpPtrVector().push_back(
1473 new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
1474 feFaceRhs->getOpPtrVector().push_back(
1475 new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
1476
1477 // create method for setting history for fluxes on boundary
1478 scaleMethodFlux = boost::shared_ptr<MethodForForceScaling>(
1479 new TimeForceScale("-flux_history", false));
1480
1481 // create method for setting history for presure heads on boundary
1482 scaleMethodValue = boost::shared_ptr<MethodForForceScaling>(
1483 new TimeForceScale("-head_history", false));
1484
1485 // Set operator to calculate essential boundary conditions
1486 feFaceBc->getOpPtrVector().push_back(
1487 new OpEvaluateBcOnFluxes(*this, "FLUXES"));
1488
1489 // Set operator to calculate initial capillary pressure
1490 feVolInitialPc->getOpPtrVector().push_back(
1491 new OpEvaluateInitiallHead(*this, "VALUES"));
1492
1493 // set residual face from Neumann terms, i.e. applied pressure
1494 feFaceRhs->getOpPtrVector().push_back(
1495 new OpRhsBcOnValues(*this, "FLUXES", scaleMethodValue));
1496 // set residual finite element operators
1497 headRateAtGaussPts = boost::make_shared<VectorDouble>();
1498 // resAtGaussPts = boost::make_shared<VectorDouble>();
1499 feVolRhs->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
1500 string("VALUES") + "_t", headRateAtGaussPts, MBTET));
1501 feVolRhs->getOpPtrVector().push_back(
1502 new OpValuesAtGaussPts(*this, "VALUES"));
1503 feVolRhs->getOpPtrVector().push_back(
1504 new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
1505 feVolRhs->getOpPtrVector().push_back(new OpResidualFlux(*this, "FLUXES"));
1506 feVolRhs->getOpPtrVector().push_back(new OpResidualMass(*this, "VALUES"));
1507 feVolRhs->getOpPtrVector().back().opType =
1508 ForcesAndSourcesCore::UserDataOperator::OPROW;
1509 // set tangent matrix finite element operators
1510 feVolLhs->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
1511 string("VALUES") + "_t", headRateAtGaussPts, MBTET));
1512 // feVolLhs->getOpPtrVector().push_back(
1513 // new
1514 // OpCalculateScalarFieldValues(string("FLUXES")+"_residual",resAtGaussPts,MBTET)
1515 // );
1516 feVolLhs->getOpPtrVector().push_back(
1517 new OpValuesAtGaussPts(*this, "VALUES"));
1518 feVolLhs->getOpPtrVector().push_back(
1519 new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
1520 feVolLhs->getOpPtrVector().push_back(
1521 new OpTauDotSigma_HdivHdiv(*this, "FLUXES"));
1522 feVolLhs->getOpPtrVector().push_back(new OpVU_L2L2(*this, "VALUES"));
1523 feVolLhs->getOpPtrVector().push_back(
1524 new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES"));
1525 feVolLhs->getOpPtrVector().push_back(
1526 new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES"));
1527
1528 // Adding finite elements to DM, time solver will ask for it to assemble
1529 // tangent matrices and residuals
1530 boost::shared_ptr<FEMethod> null;
1531 CHKERR DMMoFEMTSSetIFunction(dM, "MIX_BCVALUE", feFaceRhs, null, null);
1532 CHKERR DMMoFEMTSSetIFunction(dM, "MIX", feVolRhs, null, null);
1533 CHKERR DMMoFEMTSSetIJacobian(dM, "MIX", feVolLhs, null, null);
1534
1535 // setting up post-processing
1536
1537
1538 auto get_post_process_ele = [&]() {
1539 auto post_process = boost::make_shared<PostProcEle>(mField);
1540
1541 auto jac_ptr = boost::make_shared<MatrixDouble>();
1542 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
1543 auto det_ptr = boost::make_shared<VectorDouble>();
1544 post_process->getOpPtrVector().push_back(new OpCalculateHOJac<3>(jac_ptr));
1545 post_process->getOpPtrVector().push_back(
1546 new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
1547 post_process->getOpPtrVector().push_back(
1548 new OpSetHOContravariantPiolaTransform(HDIV, det_ptr, jac_ptr));
1549 post_process->getOpPtrVector().push_back(
1550 new OpSetHOInvJacToScalarBases<3>(L2, inv_jac_ptr));
1551
1552 auto values_ptr = boost::make_shared<VectorDouble>();
1553 auto grad_ptr = boost::make_shared<MatrixDouble>();
1554 auto flux_ptr = boost::make_shared<MatrixDouble>();
1555
1556 post_process->getOpPtrVector().push_back(
1557 new OpCalculateScalarFieldValues("VALUES", values_ptr));
1558 post_process->getOpPtrVector().push_back(
1559 new OpCalculateScalarFieldGradient<3>("VALUES", grad_ptr));
1560 post_process->getOpPtrVector().push_back(
1561 new OpCalculateHVecVectorField<3>("FLUXES", flux_ptr));
1562
1563 using OpPPMap = OpPostProcMapInMoab<3, 3>;
1564
1565 post_process->getOpPtrVector().push_back(
1566
1567 new OpPPMap(post_process->getPostProcMesh(),
1568 post_process->getMapGaussPts(),
1569
1570 {{"VALUES", values_ptr}},
1571
1572 {{"GRAD", grad_ptr}, {"FLUXES", flux_ptr}},
1573
1574 {}, {}
1575
1576 ));
1577
1578 return post_process;
1579 };
1580
1581 auto post_process = get_post_process_ele();
1582
1583 post_process->getOpPtrVector().push_back(
1584 new OpValuesAtGaussPts(*this, "VALUES"));
1585 post_process->getOpPtrVector().push_back(
1586 new OpPostProcMaterial(*this, post_process->getPostProcMesh(),
1587 post_process->getMapGaussPts(), "VALUES"));
1588
1589 // Integrate fluxes on boundary
1590 boost::shared_ptr<ForcesAndSourcesCore> flux_integrate;
1591 flux_integrate = boost::shared_ptr<ForcesAndSourcesCore>(
1593 flux_integrate->getOpPtrVector().push_back(
1594 new OpIntegrateFluxes(*this, "FLUXES"));
1595 int frequency = 1;
1596 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Monitor post proc", "none");
1597 CHKERR PetscOptionsInt(
1598 "-how_often_output",
1599 "frequency how often results are dumped on hard disk", "", frequency,
1600 &frequency, NULL);
1601 ierr = PetscOptionsEnd();
1602 CHKERRG(ierr);
1603
1604 tsMonitor = boost::shared_ptr<FEMethod>(
1605 new MonitorPostProc(*this, post_process, flux_integrate, frequency));
1606 TsCtx *ts_ctx;
1607 DMMoFEMGetTsCtx(dM, &ts_ctx);
1608 ts_ctx->getPostProcessMonitor().push_back(tsMonitor);
1610 }
1611
1612 Vec D1; ///< Vector with inital head capillary pressure
1613 Vec ghostFlux; ///< Ghost Vector of integrated fluxes
1614
1615 /**
1616 * \brief Create vectors and matrices
1617 * @return Error code
1618 */
1621 CHKERR DMCreateMatrix(dM, &Aij);
1622 CHKERR DMCreateGlobalVector(dM, &D0);
1623 CHKERR VecDuplicate(D0, &D1);
1624 CHKERR VecDuplicate(D0, &D);
1625 CHKERR VecDuplicate(D0, &F);
1626 int ghosts[] = {0};
1627 int nb_locals = mField.get_comm_rank() == 0 ? 1 : 0;
1628 int nb_ghosts = mField.get_comm_rank() > 0 ? 1 : 0;
1629 CHKERR VecCreateGhost(PETSC_COMM_WORLD, nb_locals, 1, nb_ghosts, ghosts,
1630 &ghostFlux);
1632 }
1633
1634 /**
1635 * \brief Delete matrices and vector when no longer needed
1636 * @return error code
1637 */
1640 CHKERR MatDestroy(&Aij);
1641 CHKERR VecDestroy(&D);
1642 CHKERR VecDestroy(&D0);
1643 CHKERR VecDestroy(&D1);
1644 CHKERR VecDestroy(&F);
1645 CHKERR VecDestroy(&ghostFlux);
1647 }
1648
1649 /**
1650 * \brief Calculate boundary conditions for fluxes
1651 * @return Error code
1652 */
1655 // clear vectors
1656 CHKERR VecZeroEntries(D0);
1657 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
1658 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
1659 // clear essential bc indices, it could have dofs from other mesh refinement
1660 bcIndices.clear();
1661 // set operator to calculate essential boundary conditions
1662 CHKERR DMoFEMLoopFiniteElements(dM, "MIX_BCFLUX", feFaceBc);
1663 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_REVERSE);
1664 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_REVERSE);
1665 CHKERR VecAssemblyBegin(D0);
1666 CHKERR VecAssemblyEnd(D0);
1667 double norm2D0;
1668 CHKERR VecNorm(D0, NORM_2, &norm2D0);
1669 // CHKERR VecView(D0,PETSC_VIEWER_STDOUT_WORLD);
1670 PetscPrintf(PETSC_COMM_WORLD, "norm2D0 = %6.4e\n", norm2D0);
1672 }
1673
1674 /**
1675 * \brief Calculate inital pressure head distribution
1676 * @return Error code
1677 */
1680 // clear vectors
1681 CHKERR VecZeroEntries(D1);
1682 CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_FORWARD);
1683 CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_FORWARD);
1684 // Calculate initial pressure head on each element
1685 CHKERR DMoFEMLoopFiniteElements(dM, "MIX", feVolInitialPc);
1686 // Assemble vector
1687 CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_REVERSE);
1688 CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_REVERSE);
1689 CHKERR VecAssemblyBegin(D1);
1690 CHKERR VecAssemblyEnd(D1);
1691 // Calculate norm
1692 double norm2D1;
1693 CHKERR VecNorm(D1, NORM_2, &norm2D1);
1694 // CHKERR VecView(D0,PETSC_VIEWER_STDOUT_WORLD);
1695 PetscPrintf(PETSC_COMM_WORLD, "norm2D1 = %6.4e\n", norm2D1);
1697 }
1698
1699 /**
1700 * \brief solve problem
1701 * @return error code
1702 */
1703 MoFEMErrorCode solveProblem(bool set_initial_pc = true) {
1705 if (set_initial_pc) {
1706 // Set initial head
1707 CHKERR DMoFEMMeshToLocalVector(dM, D1, INSERT_VALUES, SCATTER_REVERSE);
1708 }
1709
1710 // Initiate vector from data on the mesh
1711 CHKERR DMoFEMMeshToLocalVector(dM, D, INSERT_VALUES, SCATTER_FORWARD);
1712
1713 // Create time solver
1714 TS ts;
1715 CHKERR TSCreate(PETSC_COMM_WORLD, &ts);
1716 // Use backward Euler method
1717 CHKERR TSSetType(ts, TSBEULER);
1718 // Set final time
1719 double ftime = 1;
1720 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
1721 // Setup solver from commabd line
1722 CHKERR TSSetFromOptions(ts);
1723 // Set DM to TS
1724 CHKERR TSSetDM(ts, dM);
1725#if PETSC_VERSION_GE(3, 7, 0)
1726 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
1727#endif
1728 // Set-up monitor
1729 TsCtx *ts_ctx;
1730 DMMoFEMGetTsCtx(dM, &ts_ctx);
1731 CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx, PETSC_NULL);
1732
1733 // This add SNES monitor, to show error by fields. It is dirty trick
1734 // to add monitor, so code is hidden from doxygen
1735 CHKERR TSSetSolution(ts, D);
1736 CHKERR TSSetUp(ts);
1737 SNES snes;
1738 CHKERR TSGetSNES(ts, &snes);
1739
1740#if PETSC_VERSION_GE(3, 7, 0)
1741 {
1742 PetscViewerAndFormat *vf;
1743 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
1744 PETSC_VIEWER_DEFAULT, &vf);
1745 CHKERR SNESMonitorSet(
1746 snes,
1747 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
1748 void *))SNESMonitorFields,
1749 vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
1750 }
1751#else
1752 {
1753 CHKERR SNESMonitorSet(snes,
1754 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
1755 void *))SNESMonitorFields,
1756 0, 0);
1757 }
1758#endif
1759
1760 CHKERR TSSolve(ts, D);
1761
1762 // Get statistic form TS and print it
1763 CHKERR TSGetTime(ts, &ftime);
1764 PetscInt steps, snesfails, rejects, nonlinits, linits;
1765 CHKERR TSGetTimeStepNumber(ts, &steps);
1766 CHKERR TSGetSNESFailures(ts, &snesfails);
1767 CHKERR TSGetStepRejections(ts, &rejects);
1768 CHKERR TSGetSNESIterations(ts, &nonlinits);
1769 CHKERR TSGetKSPIterations(ts, &linits);
1770 PetscPrintf(PETSC_COMM_WORLD,
1771 "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits "
1772 "%D, linits %D\n",
1773 steps, rejects, snesfails, ftime, nonlinits, linits);
1774
1776 }
1777};
1778
1779} // namespace MixTransport
1780
1781#endif // __UNSATURATD_FLOW_HPP__
static Index< 'K', 3 > K
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
EntitiesFieldData::EntData EntData
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
Definition: cholesky.hpp:221
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
Definition: cholesky.hpp:52
@ ROW
Definition: definitions.h:123
@ MF_ZERO
Definition: definitions.h:98
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:1070
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
Definition: DMMMoFEM.cpp:747
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: DMMMoFEM.cpp:118
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:450
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:470
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMMoFEM.cpp:1089
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:533
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, 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 TS Jacobian evaluation function
Definition: DMMMoFEM.cpp:800
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
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
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
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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.
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.
#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.
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
const FTensor::Tensor2< T, Dim, Dim > Vec
auto f
Definition: HenckyOps.hpp:5
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:217
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
double w(const double g, const double t)
const double D
diffusivity
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double scale
Definition: plastic.cpp:94
constexpr auto field_name
Generic material model for unsaturated water transport.
virtual MoFEMErrorCode calC()
virtual MoFEMErrorCode calSe()
virtual void printMatParameters(const int id, const std::string &prefix) const =0
virtual MoFEMErrorCode calTheta()
double diffC
Derivative of capacity [S^2/L^2 * L^2/F ].
static double ePsilon0
Regularization parameter.
virtual MoFEMErrorCode calK()
virtual MoFEMErrorCode calDiffC()
double sCale
Scale time dependent eq.
double K
Hydraulic conductivity [L/s].
Range tEts
Elements with this material.
double C
Capacity [S^2/L^2].
double Se
Effective saturation.
static double scaleZ
Scale z direction.
virtual double initialPcEval() const =0
Initialize head.
double h_t
rate of hydraulic head
static double ePsilon1
Regularization parameter.
virtual MoFEMErrorCode calDiffK()
double diffK
Derivative of hydraulic conductivity [L/s * L^2/F].
VectorDouble valuesAtGaussPts
values at integration points on element
VectorDouble divergenceAtGaussPts
divergence at integration points on element
MatrixDouble fluxesAtGaussPts
fluxes at integration points on element
Class storing information about boundary condition.
boost::function< double(const double x, const double y, const double z)> hookFun
Set integration rule to boundary elements.
int operator()(int p_row, int p_col, int p_data) const
MoFEMErrorCode operator()()
function is run for every finite element
MonitorPostProc(UnsaturatedFlowElement &ctx, boost::shared_ptr< PostProcEle > &post_proc, boost::shared_ptr< ForcesAndSourcesCore > flux_Integrate, const int frequency)
boost::shared_ptr< PostProcEle > postProc
boost::shared_ptr< ForcesAndSourcesCore > fluxIntegrate
MoFEMErrorCode preProcess()
function is run at the beginning of loop
UnsaturatedFlowElement & cTx
MoFEMErrorCode postProcess()
function is run at the end of loop
const int fRequency
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations.
OpDivTauU_HdivL2(UnsaturatedFlowElement &ctx, const std::string &flux_name_col, const std::string &val_name_row)
Constructor.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpEvaluateInitiallHead(UnsaturatedFlowElement &ctx, const std::string &val_name)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrate boundary condition.
OpIntegrateFluxes(UnsaturatedFlowElement &ctx, const std::string fluxes_name)
Constructor.
OpPostProcMaterial(UnsaturatedFlowElement &ctx, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, const std::string field_name)
moab::Interface & postProcMesh
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
UnsaturatedFlowElement & cTx
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpResidualFlux(UnsaturatedFlowElement &ctx, const std::string &flux_name)
OpResidualMass(UnsaturatedFlowElement &ctx, const std::string &val_name)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Evaluate boundary condition at the boundary.
OpRhsBcOnValues(UnsaturatedFlowElement &ctx, const std::string fluxes_name, boost::shared_ptr< MethodForForceScaling > &value_scale)
Constructor.
boost::shared_ptr< MethodForForceScaling > valueScale
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrate boundary condition.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble matrix.
OpTauDotSigma_HdivHdiv(UnsaturatedFlowElement &ctx, const std::string flux_name)
OpVDivSigma_L2Hdiv(UnsaturatedFlowElement &ctx, const std::string &val_name_row, const std::string &flux_name_col)
Constructor.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations.
OpVU_L2L2(UnsaturatedFlowElement &ctx, const std::string value_name)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble matrix.
Set integration rule to volume elements.
int operator()(int, int, int p_data) const
Post proces method for volume element Assemble vectors and matrices and apply essential boundary cond...
postProcessVol(UnsaturatedFlowElement &ctx, boost::shared_ptr< ForcesAndSourcesCore > &fe_ptr)
boost::shared_ptr< ForcesAndSourcesCore > fePtr
Pre-peprocessing Set head pressute rate and get inital essential boundary conditions.
boost::shared_ptr< ForcesAndSourcesCore > fePtr
preProcessVol(UnsaturatedFlowElement &ctx, boost::shared_ptr< ForcesAndSourcesCore > &fe_ptr)
Implementation of operators, problem and finite elements for unsaturated flow.
MoFEMErrorCode calculateEssentialBc()
Calculate boundary conditions for fluxes.
boost::shared_ptr< ForcesAndSourcesCore > feVolLhs
Assemble tangent matrix.
MoFEMErrorCode getBcOnValues(const EntityHandle ent, const int gg, const double x, const double y, const double z, double &value)
Get value on boundary.
std::map< int, boost::shared_ptr< GenericMaterial > > MaterialsDoubleMap
MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x, const double y, const double z, double &flux)
essential (Neumann) boundary condition (set fluxes)
MoFEMErrorCode setFiniteElements(ForcesAndSourcesCore::RuleHookFun vol_rule=VolRule(), ForcesAndSourcesCore::RuleHookFun face_rule=FaceRule())
Create finite element instances.
boost::shared_ptr< VectorDouble > headRateAtGaussPts
Vector keeps head rate.
boost::shared_ptr< MethodForForceScaling > scaleMethodFlux
Method scaling fluxes.
MoFEMErrorCode buildProblem(Range zero_flux_ents, BitRefLevel ref_level=BitRefLevel().set(0))
Build problem.
Vec D1
Vector with inital head capillary pressure.
map< int, boost::shared_ptr< BcData > > BcMap
MaterialsDoubleMap dMatMap
materials database
boost::shared_ptr< ForcesAndSourcesCore > feFaceRhs
Face element apply natural bc.
Vec ghostFlux
Ghost Vector of integrated fluxes.
MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes, const int order)
add fields
boost::shared_ptr< ForcesAndSourcesCore > feVolRhs
Assemble residual vector.
UnsaturatedFlowElement(MoFEM::Interface &m_field)
boost::shared_ptr< FEMethod > tsMonitor
BcMap bcValueMap
Store boundary condition for head capillary pressure.
MoFEMErrorCode solveProblem(bool set_initial_pc=true)
solve problem
MoFEMErrorCode destroyMatrices()
Delete matrices and vector when no longer needed.
boost::shared_ptr< ForcesAndSourcesCore > feVolInitialPc
Calculate inital boundary conditions.
MoFEMErrorCode calculateInitialPc()
Calculate inital pressure head distribution.
DM dM
Discrete manager for unsaturated flow problem.
MoFEMErrorCode createMatrices()
Create vectors and matrices.
boost::shared_ptr< ForcesAndSourcesCore > feFaceBc
Elemnet to calculate essential bc.
MoFEMErrorCode addFiniteElements(const std::string &fluxes_name, const std::string &values_name)
add finite elements
virtual MoFEMErrorCode getMaterial(const EntityHandle ent, int &block_id) const
For given element handle get material block Id.
boost::shared_ptr< MethodForForceScaling > scaleMethodValue
Method scaling values.
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MPI_Comm & get_comm() const =0
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.
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
structure for User Loop Methods on finite elements
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
auto getFTensor0IntegrationWeight()
Get integration weights.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Post post-proc data at points from hash maps.
TS ts
time solver
PetscReal ts_t
time
Vec & ts_F
residual vector
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
Mat & ts_B
Preconditioner for ts_A.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Force scale operator for reading two columns.