v0.13.2
Loading...
Searching...
No Matches
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
61 virtual MoFEMErrorCode calDiffK() {
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
76 virtual MoFEMErrorCode calDiffC() {
78 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
79 "Not implemented how to calculate capacity");
81 }
82
83 virtual MoFEMErrorCode calTheta() {
85 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
86 "Not implemented how to calculate capacity");
88 }
89
90 virtual MoFEMErrorCode calSe() {
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
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
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 */
265 MoFEMErrorCode doWork(int side, EntityType type,
266 EntitiesFieldData::EntData &data) {
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
314 VectorDouble divVec, nF;
316
317 MoFEMErrorCode doWork(int side, EntityType type,
318 EntitiesFieldData::EntData &data) {
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
336 auto t_h = getFTensor0FromVec(cTx.valuesAtGaussPts);
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
403 VectorDouble nF;
404
405 MoFEMErrorCode doWork(int side, EntityType type,
406 EntitiesFieldData::EntData &data) {
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
422 auto t_h = getFTensor0FromVec(cTx.valuesAtGaussPts);
423 // Get pressure rate
424 auto t_h_t = getFTensor0FromVec(*cTx.headRateAtGaussPts);
425 // Flux divergence
426 auto t_div_flux = getFTensor0FromVec(cTx.divergenceAtGaussPts);
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
474 MatrixDouble nN, transNN;
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,
490 EntitiesFieldData::EntData &row_data,
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
509 auto t_h = getFTensor0FromVec(cTx.valuesAtGaussPts);
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
571 MatrixDouble nN;
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,
585 EntitiesFieldData::EntData &row_data,
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
604 auto t_h = getFTensor0FromVec(cTx.valuesAtGaussPts);
605 // // Get pressure
606 // auto t_flux_residual = getFTensor0FromVec(*cTx.resAtGaussPts);
607 // Get pressure rate
608 auto t_h_t = getFTensor0FromVec(*cTx.headRateAtGaussPts);
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
681 MatrixDouble nN;
682 VectorDouble divVec;
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,
698 EntitiesFieldData::EntData &row_data,
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
753 MatrixDouble nN;
754 VectorDouble divVec;
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,
771 EntitiesFieldData::EntData &row_data,
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
791 auto t_h = getFTensor0FromVec(cTx.valuesAtGaussPts);
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
855 MatrixDouble nN;
856 VectorDouble nF;
857
858 MoFEMErrorCode doWork(int side, EntityType type,
859 EntitiesFieldData::EntData &data) {
861 if (data.getFieldData().size() == 0)
863 auto nb_dofs = data.getFieldData().size();
864 if (nb_dofs != data.getN().size2()) {
865 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
866 "wrong number of dofs");
867 }
868 nN.resize(nb_dofs, nb_dofs, false);
869 nF.resize(nb_dofs, false);
870 nN.clear();
871 nF.clear();
872
873 // Get EntityHandle of the finite element
875
876
877 // Get material block id
878 int block_id;
879 CHKERR cTx.getMaterial(fe_ent, block_id);
880 // Get material block
881 boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
882
883 // loop over integration points
884 auto nb_gauss_pts = getGaussPts().size2();
885 for (auto gg = 0; gg != nb_gauss_pts; gg++) {
886 // get coordinates at integration point
887 block_data->x = getCoordsAtGaussPts()(gg, 0);
888 block_data->y = getCoordsAtGaussPts()(gg, 1);
889 block_data->z = getCoordsAtGaussPts()(gg, 2);
890 // get weight for integration rule
891 double alpha = getGaussPts()(2, gg);
892 nN += alpha * outer_prod(data.getN(gg), data.getN(gg));
893 nF += alpha * block_data->initialPcEval() * data.getN(gg);
894 }
895
896 // factor matrix
898 // solve local problem
899 cholesky_solve(nN, nF, ublas::lower());
900
901 // set solution to vector
902 CHKERR VecSetValues(cTx.D1, nb_dofs, &*data.getIndices().begin(),
903 &*nF.begin(), INSERT_VALUES);
904
906 }
907 };
908
912
913 /**
914 * \brief Constructor
915 */
917 const std::string fluxes_name)
919 fluxes_name, UserDataOperator::OPROW),
920 cTx(ctx) {}
921
923
924 /**
925 * \brief Integrate boundary condition
926 * @param side local index of entity
927 * @param type type of entity
928 * @param data data on entity
929 * @return error code
930 */
931 MoFEMErrorCode doWork(int side, EntityType type,
932 EntitiesFieldData::EntData &data) {
934 int nb_dofs = data.getFieldData().size();
935 if (nb_dofs == 0)
937 // Get base function
938 auto t_n_hdiv = data.getFTensor1N<3>();
939 // get normal of face
940 auto t_normal = getFTensor1Normal();
941 // Integration weight
942 auto t_w = getFTensor0IntegrationWeight();
943 double flux_on_entity = 0;
944 int nb_gauss_pts = data.getN().size1();
945 for (int gg = 0; gg < nb_gauss_pts; gg++) {
946 auto t_data = data.getFTensor0FieldData();
947 for (int rr = 0; rr != nb_dofs; rr++) {
948 flux_on_entity -= (0.5 * t_data * t_w) * (t_n_hdiv(i) * t_normal(i));
949 ++t_n_hdiv;
950 ++t_data;
951 }
952 ++t_w;
953 }
954 CHKERR VecSetValue(cTx.ghostFlux, 0, flux_on_entity, ADD_VALUES);
956 }
957 };
958
959 /**
960 * Operator used to post-process results for unsaturated infiltration problem.
961 * Operator should with element for post-processing results, i.e.
962 * PostProcVolumeOnRefinedMesh
963 */
966
968 moab::Interface &postProcMesh;
969 std::vector<EntityHandle> &mapGaussPts;
970
972 moab::Interface &post_proc_mesh,
973 std::vector<EntityHandle> &map_gauss_pts,
974 const std::string field_name)
975 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
977 cTx(ctx), postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts) {}
978
979 MoFEMErrorCode doWork(int side, EntityType type,
980 EntitiesFieldData::EntData &data) {
982 int nb_dofs = data.getFieldData().size();
983 if (nb_dofs == 0)
985
986 // Get EntityHandle of the finite element
987 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
988 // Get material block id
989 int block_id;
990 CHKERR cTx.getMaterial(fe_ent, block_id);
991 // Get material block
992 boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
993
994 // Set bloc Id
995 Tag th_id;
996 int def_block_id = -1;
997 CHKERR postProcMesh.tag_get_handle("BLOCK_ID", 1, MB_TYPE_INTEGER, th_id,
998 MB_TAG_CREAT | MB_TAG_SPARSE,
999 &def_block_id);
1000
1001 // Create mesh tag. Tags are created on post-processing mesh and
1002 // visable in post-processor, e.g. Paraview
1003 double zero = 0;
1004 Tag th_theta;
1005 CHKERR postProcMesh.tag_get_handle("THETA", 1, MB_TYPE_DOUBLE, th_theta,
1006 MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
1007 Tag th_se;
1008 CHKERR postProcMesh.tag_get_handle("Se", 1, MB_TYPE_DOUBLE, th_se,
1009 MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
1010 // Tag th_ks;
1011 // CHKERR postProcMesh.tag_get_handle(
1012 // "Ks",1,MB_TYPE_DOUBLE,th_ks,
1013 // MB_TAG_CREAT|MB_TAG_SPARSE,&zero
1014 // );
1015 // CHKERR postProcMesh.tag_set_data(th_ks,&fe_ent,1,&block_data->Ks);
1016 Tag th_k;
1017 CHKERR postProcMesh.tag_get_handle("K", 1, MB_TYPE_DOUBLE, th_k,
1018 MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
1019 Tag th_c;
1020 CHKERR postProcMesh.tag_get_handle("C", 1, MB_TYPE_DOUBLE, th_c,
1021 MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
1022
1023 // Get pressure at integration points
1024 auto t_h = getFTensor0FromVec(cTx.valuesAtGaussPts);
1025 // Coords at integration points
1026 auto t_coords = getFTensor1CoordsAtGaussPts();
1027
1028 int nb_gauss_pts = data.getN().size1();
1029 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1030 block_data->h = t_h;
1031 block_data->x = t_coords(0);
1032 block_data->y = t_coords(1);
1033 block_data->z = t_coords(2);
1034 // Calculate theta (water content) and save it on mesh tags
1035 CHKERR block_data->calTheta();
1036 double theta = block_data->tHeta;
1037 CHKERR postProcMesh.tag_set_data(th_theta, &mapGaussPts[gg], 1, &theta);
1038 CHKERR block_data->calSe();
1039 // Calculate Se (effective saturation and save it on the mesh tags)
1040 double Se = block_data->Se;
1041 CHKERR postProcMesh.tag_set_data(th_se, &mapGaussPts[gg], 1, &Se);
1042 // Calculate K (hydraulic conductivity) and save it on the mesh tags
1043 CHKERR block_data->calK();
1044 double K = block_data->K;
1045 CHKERR postProcMesh.tag_set_data(th_k, &mapGaussPts[gg], 1, &K);
1046 // Calculate water capacity and save it on the mesh tags
1047 CHKERR block_data->calC();
1048 double C = block_data->C;
1049 CHKERR postProcMesh.tag_set_data(th_c, &mapGaussPts[gg], 1, &C);
1050
1051 // Set block iD
1052 CHKERR postProcMesh.tag_set_data(th_id, &mapGaussPts[gg], 1, &block_id);
1053
1054 ++t_h;
1055 ++t_coords;
1056 }
1057
1059 }
1060 };
1061
1062 /**
1063 * Finite element implementation called by TS monitor. Element calls other
1064 * finite elements to evaluate material properties and save results on the
1065 * mesh.
1066 *
1067 * \note Element overloaded only FEMethod::postProcess methods where other
1068 * elements are called.
1069 */
1071
1073 boost::shared_ptr<PostProcEle> postProc;
1074 boost::shared_ptr<ForcesAndSourcesCore> fluxIntegrate;
1075
1076 const int fRequency;
1077
1079 boost::shared_ptr<PostProcEle> &post_proc,
1080 boost::shared_ptr<ForcesAndSourcesCore> flux_Integrate,
1081 const int frequency)
1082 : cTx(ctx), postProc(post_proc), fluxIntegrate(flux_Integrate),
1083 fRequency(frequency) {}
1084
1085 MoFEMErrorCode preProcess() {
1088 }
1089
1090 MoFEMErrorCode operator()() {
1093 }
1094
1095 MoFEMErrorCode postProcess() {
1097
1098 // Get time step
1099 int step;
1100 CHKERR TSGetTimeStepNumber(ts, &step);
1101
1102 if ((step) % fRequency == 0) {
1103 // Post-process results and save in the file
1104 PetscPrintf(PETSC_COMM_WORLD, "Output results %d - %d\n", step,
1105 fRequency);
1106 CHKERR DMoFEMLoopFiniteElements(cTx.dM, "MIX", postProc);
1107 CHKERR postProc->writeFile(
1108 string("out_") + boost::lexical_cast<std::string>(step) + ".h5m");
1109 }
1110
1111 // Integrate fluxes on faces where pressure head is applied
1112 CHKERR VecZeroEntries(cTx.ghostFlux);
1113 CHKERR VecGhostUpdateBegin(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
1114 CHKERR VecGhostUpdateEnd(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
1115 // Run finite element to integrate fluxes
1116 CHKERR DMoFEMLoopFiniteElements(cTx.dM, "MIX_BCVALUE", fluxIntegrate);
1117 CHKERR VecAssemblyBegin(cTx.ghostFlux);
1118 CHKERR VecAssemblyEnd(cTx.ghostFlux);
1119 // accumulate errors from processors
1120 CHKERR VecGhostUpdateBegin(cTx.ghostFlux, ADD_VALUES, SCATTER_REVERSE);
1121 CHKERR VecGhostUpdateEnd(cTx.ghostFlux, ADD_VALUES, SCATTER_REVERSE);
1122 // scatter errors to all processors
1123 CHKERR VecGhostUpdateBegin(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
1124 CHKERR VecGhostUpdateEnd(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
1125 double *ghost_flux;
1126 CHKERR VecGetArray(cTx.ghostFlux, &ghost_flux);
1127 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Flux at time %6.4g %6.4g\n", ts_t,
1128 ghost_flux[0]);
1129 CHKERR VecRestoreArray(cTx.ghostFlux, &ghost_flux);
1130
1132 }
1133 };
1134
1135 /// \brief add fields
1136 MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes,
1137 const int order) {
1139 // Fields
1142 CHKERR mField.add_field(values + "_t", L2, AINSWORTH_LEGENDRE_BASE, 1);
1143 // CHKERR mField.add_field(fluxes+"_residual",L2,AINSWORTH_LEGENDRE_BASE,1);
1144
1145 // meshset consisting all entities in mesh
1146 EntityHandle root_set = mField.get_moab().get_root_set();
1147 // add entities to field
1148
1150 if (it->getName().compare(0, 4, "SOIL") != 0)
1151 continue;
1152 CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
1153 MBTET, fluxes);
1154 CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
1155 MBTET, values);
1156 CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
1157 MBTET, values + "_t");
1158 // CHKERR mField.add_ents_to_field_by_type(
1159 // dMatMap[it->getMeshsetId()]->tEts,MBTET,fluxes+"_residual"
1160 // );
1161 }
1162
1163 CHKERR mField.set_field_order(root_set, MBTET, fluxes, order + 1);
1164 CHKERR mField.set_field_order(root_set, MBTRI, fluxes, order + 1);
1165 CHKERR mField.set_field_order(root_set, MBTET, values, order);
1166 CHKERR mField.set_field_order(root_set, MBTET, values + "_t", order);
1167 // CHKERR mField.set_field_order(root_set,MBTET,fluxes+"_residual",order);
1169 }
1170
1171 /// \brief add finite elements
1172 MoFEMErrorCode addFiniteElements(const std::string &fluxes_name,
1173 const std::string &values_name) {
1175
1176 // Define element "MIX". Note that this element will work with fluxes_name
1177 // and values_name. This reflect bilinear form for the problem
1186 values_name + "_t");
1187 // CHKERR
1188 // mField.modify_finite_element_add_field_data("MIX",fluxes_name+"_residual");
1189
1191 if (it->getName().compare(0, 4, "SOIL") != 0)
1192 continue;
1194 dMatMap[it->getMeshsetId()]->tEts, MBTET, "MIX");
1195 }
1196
1197 // Define element to integrate natural boundary conditions, i.e. set values.
1198 CHKERR mField.add_finite_element("MIX_BCVALUE", MF_ZERO);
1200 fluxes_name);
1202 fluxes_name);
1204 fluxes_name);
1206 values_name);
1207
1209 if (it->getName().compare(0, 4, "HEAD") != 0)
1210 continue;
1212 bcValueMap[it->getMeshsetId()]->eNts, MBTRI, "MIX_BCVALUE");
1213 }
1214
1215 // Define element to apply essential boundary conditions.
1216 CHKERR mField.add_finite_element("MIX_BCFLUX", MF_ZERO);
1218 fluxes_name);
1220 fluxes_name);
1222 fluxes_name);
1224 values_name);
1225
1227 if (it->getName().compare(0, 4, "FLUX") != 0)
1228 continue;
1230 bcFluxMap[it->getMeshsetId()]->eNts, MBTRI, "MIX_BCFLUX");
1231 }
1232
1234 }
1235
1236 /**
1237 * \brief Build problem
1238 * @param ref_level mesh refinement on which mesh problem you like to built.
1239 * @return error code
1240 */
1241 MoFEMErrorCode buildProblem(Range zero_flux_ents,
1242 BitRefLevel ref_level = BitRefLevel().set(0)) {
1244
1245 // Build fields
1247 // Build finite elements
1249 CHKERR mField.build_finite_elements("MIX_BCFLUX");
1250 CHKERR mField.build_finite_elements("MIX_BCVALUE");
1251 // Build adjacencies of degrees of freedom and elements
1252 CHKERR mField.build_adjacencies(ref_level);
1253
1254 // create DM instance
1255 CHKERR DMCreate(PETSC_COMM_WORLD, &dM);
1256 // setting that DM is type of DMMOFEM, i.e. MOFEM implementation manages DM
1257 CHKERR DMSetType(dM, "DMMOFEM");
1258 // mesh is portioned, each process keeps only part of problem
1259 CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
1260 // creates problem in DM
1261 CHKERR DMMoFEMCreateMoFEM(dM, &mField, "MIX", ref_level);
1262 // discretised problem creates square matrix (that makes some optimizations)
1263 CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
1264 // set DM options from command line
1265 CHKERR DMSetFromOptions(dM);
1266 // add finite elements
1267 CHKERR DMMoFEMAddElement(dM, "MIX");
1268 CHKERR DMMoFEMAddElement(dM, "MIX_BCFLUX");
1269 CHKERR DMMoFEMAddElement(dM, "MIX_BCVALUE");
1270 // constructor data structures
1271 CHKERR DMSetUp(dM);
1272
1273 // remove zero flux dofs
1274 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
1275 "MIX", "FLUXES", zero_flux_ents);
1276
1277 PetscSection section;
1278 CHKERR mField.getInterface<ISManager>()->sectionCreate("MIX", &section);
1279 CHKERR DMSetSection(dM, section);
1280 // CHKERR PetscSectionView(section,PETSC_VIEWER_STDOUT_WORLD);
1281 CHKERR PetscSectionDestroy(&section);
1282
1284 }
1285
1286 boost::shared_ptr<ForcesAndSourcesCore>
1287 feFaceBc; ///< Elemnet to calculate essential bc
1288 boost::shared_ptr<ForcesAndSourcesCore>
1289 feFaceRhs; ///< Face element apply natural bc
1290 boost::shared_ptr<ForcesAndSourcesCore>
1291 feVolInitialPc; ///< Calculate inital boundary conditions
1292 boost::shared_ptr<ForcesAndSourcesCore>
1293 feVolRhs; ///< Assemble residual vector
1294 boost::shared_ptr<ForcesAndSourcesCore> feVolLhs; ///< Assemble tangent matrix
1295 boost::shared_ptr<MethodForForceScaling>
1296 scaleMethodFlux; ///< Method scaling fluxes
1297 boost::shared_ptr<MethodForForceScaling>
1298 scaleMethodValue; ///< Method scaling values
1299 boost::shared_ptr<FEMethod> tsMonitor; ///< Element used by TS monitor to
1300 ///< postprocess results at time step
1301
1302 boost::shared_ptr<VectorDouble>
1303 headRateAtGaussPts; ///< Vector keeps head rate
1304 // boost::shared_ptr<VectorDouble> resAtGaussPts; ///< Residual field
1305
1306 /**
1307 * \brief Set integration rule to volume elements
1308 *
1309 */
1310 struct VolRule {
1311 int operator()(int, int, int p_data) const { return 3 * p_data; }
1312 };
1313 /**
1314 * \brief Set integration rule to boundary elements
1315 *
1316 */
1317 struct FaceRule {
1318 int operator()(int p_row, int p_col, int p_data) const {
1319 return 2 * p_data;
1320 }
1321 };
1322
1323 std::vector<int> bcVecIds;
1324 VectorDouble bcVecVals, vecValsOnBc;
1325
1326 /**
1327 * \brief Pre-peprocessing
1328 * Set head pressure rate and get inital essential boundary conditions
1329 */
1332 boost::shared_ptr<ForcesAndSourcesCore> fePtr;
1333 // std::string mArk;
1334
1337 boost::shared_ptr<ForcesAndSourcesCore> &fe_ptr /*,std::string mark*/
1338 )
1339 : cTx(ctx), fePtr(fe_ptr) /*,mArk(mark)*/ {}
1340 MoFEMErrorCode operator()() {
1342 // Update pressure rates
1343 CHKERR fePtr->mField.getInterface<VecManager>()->setOtherLocalGhostVector(
1344 fePtr->problemPtr, "VALUES", string("VALUES") + "_t", ROW,
1345 fePtr->ts_u_t, INSERT_VALUES, SCATTER_REVERSE);
1346 switch (fePtr->ts_ctx) {
1347 case TSMethod::CTX_TSSETIFUNCTION:
1348 CHKERR VecSetOption(fePtr->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
1349 PETSC_TRUE);
1350 if (!cTx.bcIndices.empty()) {
1351 double scale;
1352 CHKERR cTx.scaleMethodFlux->getForceScale(fePtr->ts_t, scale);
1353 if (cTx.bcVecIds.size() != cTx.bcIndices.size()) {
1354 cTx.bcVecIds.insert(cTx.bcVecIds.begin(), cTx.bcIndices.begin(),
1355 cTx.bcIndices.end());
1356 cTx.bcVecVals.resize(cTx.bcVecIds.size(), false);
1357 cTx.vecValsOnBc.resize(cTx.bcVecIds.size(), false);
1358 }
1359 CHKERR VecGetValues(cTx.D0, cTx.bcVecIds.size(),
1360 &*cTx.bcVecIds.begin(), &*cTx.bcVecVals.begin());
1361 CHKERR VecGetValues(fePtr->ts_u, cTx.bcVecIds.size(),
1362 &*cTx.bcVecIds.begin(),
1363 &*cTx.vecValsOnBc.begin());
1364 cTx.bcVecVals *= scale;
1365 VectorDouble::iterator vit = cTx.bcVecVals.begin();
1366 const NumeredDofEntity *dof_ptr;
1367 for (std::vector<int>::iterator it = cTx.bcVecIds.begin();
1368 it != cTx.bcVecIds.end(); it++, vit++) {
1369 if (auto dof_ptr =
1370 fePtr->problemPtr->getColDofsByPetscGlobalDofIdx(*it)
1371 .lock())
1372 dof_ptr->getFieldData() = *vit;
1373 }
1374 } else {
1375 cTx.bcVecIds.resize(0);
1376 cTx.bcVecVals.resize(0);
1377 cTx.vecValsOnBc.resize(0);
1378 }
1379 break;
1380 default:
1381 // don nothing
1382 break;
1383 }
1385 }
1386 };
1387
1388 /**
1389 * \brief Post proces method for volume element
1390 * Assemble vectors and matrices and apply essential boundary conditions
1391 */
1394 boost::shared_ptr<ForcesAndSourcesCore> fePtr;
1395 // std::string mArk;
1398 boost::shared_ptr<ForcesAndSourcesCore> &fe_ptr //,std::string mark
1399 )
1400 : cTx(ctx), fePtr(fe_ptr) /*,mArk(mark)*/ {}
1401 MoFEMErrorCode operator()() {
1403 switch (fePtr->ts_ctx) {
1404 case TSMethod::CTX_TSSETIJACOBIAN: {
1405 CHKERR MatAssemblyBegin(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
1406 CHKERR MatAssemblyEnd(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
1407 // MatView(fePtr->ts_B,PETSC_VIEWER_DRAW_WORLD);
1408 // std::string wait;
1409 // std::cin >> wait;
1410 CHKERR MatZeroRowsColumns(fePtr->ts_B, cTx.bcVecIds.size(),
1411 &*cTx.bcVecIds.begin(), 1, PETSC_NULL,
1412 PETSC_NULL);
1413 CHKERR MatAssemblyBegin(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
1414 CHKERR MatAssemblyEnd(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
1415 // MatView(fePtr->ts_B,PETSC_VIEWER_DRAW_WORLD);
1416 // std::string wait;
1417 // std::cin >> wait;
1418 } break;
1419 case TSMethod::CTX_TSSETIFUNCTION: {
1420 CHKERR VecAssemblyBegin(fePtr->ts_F);
1421 CHKERR VecAssemblyEnd(fePtr->ts_F);
1422 if (!cTx.bcVecIds.empty()) {
1424 CHKERR VecSetValues(fePtr->ts_F, cTx.bcVecIds.size(),
1425 &*cTx.bcVecIds.begin(), &*cTx.vecValsOnBc.begin(),
1426 INSERT_VALUES);
1427 }
1428 CHKERR VecAssemblyBegin(fePtr->ts_F);
1429 CHKERR VecAssemblyEnd(fePtr->ts_F);
1430 } break;
1431 default:
1432 // don nothing
1433 break;
1434 }
1436 }
1437 };
1438
1439 /**
1440 * \brief Create finite element instances
1441 * @param vol_rule integration rule for volume element
1442 * @param face_rule integration rule for boundary element
1443 * @return error code
1444 */
1445 MoFEMErrorCode
1446 setFiniteElements(ForcesAndSourcesCore::RuleHookFun vol_rule = VolRule(),
1447 ForcesAndSourcesCore::RuleHookFun face_rule = FaceRule()) {
1449
1450 // create finite element instances
1451 feFaceBc = boost::shared_ptr<ForcesAndSourcesCore>(
1453 feFaceRhs = boost::shared_ptr<ForcesAndSourcesCore>(
1455 feVolInitialPc = boost::shared_ptr<ForcesAndSourcesCore>(
1456 new VolumeElementForcesAndSourcesCore(mField));
1457 feVolLhs = boost::shared_ptr<ForcesAndSourcesCore>(
1458 new VolumeElementForcesAndSourcesCore(mField));
1459 feVolRhs = boost::shared_ptr<ForcesAndSourcesCore>(
1460 new VolumeElementForcesAndSourcesCore(mField));
1461 // set integration rule to elements
1462 feFaceBc->getRuleHook = face_rule;
1463 feFaceRhs->getRuleHook = face_rule;
1464 feVolInitialPc->getRuleHook = vol_rule;
1465 feVolLhs->getRuleHook = vol_rule;
1466 feVolRhs->getRuleHook = vol_rule;
1467 // set function hook for finite element postprocessing stage
1468 feVolRhs->preProcessHook = preProcessVol(*this, feVolRhs);
1469 feVolLhs->preProcessHook = preProcessVol(*this, feVolLhs);
1470 feVolRhs->postProcessHook = postProcessVol(*this, feVolRhs);
1471 feVolLhs->postProcessHook = postProcessVol(*this, feVolLhs);
1472
1473 // Set Piola transform
1474 CHKERR AddHOOps<2, 3, 3>::add(feFaceBc->getOpPtrVector(), {HDIV});
1475 CHKERR AddHOOps<2, 3, 3>::add(feFaceRhs->getOpPtrVector(), {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 CHKERR AddHOOps<3, 3, 3>::add(feVolInitialPc->getOpPtrVector(), {HDIV, L2});
1491 feVolInitialPc->getOpPtrVector().push_back(
1492 new OpEvaluateInitiallHead(*this, "VALUES"));
1493
1494 // set residual face from Neumann terms, i.e. applied pressure
1495 feFaceRhs->getOpPtrVector().push_back(
1496 new OpRhsBcOnValues(*this, "FLUXES", scaleMethodValue));
1497 // set residual finite element operators
1498 headRateAtGaussPts = boost::make_shared<VectorDouble>();
1499
1500
1501 // resAtGaussPts = boost::make_shared<VectorDouble>();
1502 CHKERR AddHOOps<3, 3, 3>::add(feVolRhs->getOpPtrVector(), {HDIV, L2});
1503 feVolRhs->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
1504 string("VALUES") + "_t", headRateAtGaussPts, MBTET));
1505 feVolRhs->getOpPtrVector().push_back(
1506 new OpValuesAtGaussPts(*this, "VALUES"));
1507 feVolRhs->getOpPtrVector().push_back(
1508 new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
1509 feVolRhs->getOpPtrVector().push_back(new OpResidualFlux(*this, "FLUXES"));
1510 feVolRhs->getOpPtrVector().push_back(new OpResidualMass(*this, "VALUES"));
1511 feVolRhs->getOpPtrVector().back().opType =
1512 ForcesAndSourcesCore::UserDataOperator::OPROW;
1513 // set tangent matrix finite element operators
1514 CHKERR AddHOOps<3, 3, 3>::add(feVolLhs->getOpPtrVector(), {HDIV, L2});
1515 feVolLhs->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
1516 string("VALUES") + "_t", headRateAtGaussPts, MBTET));
1517 // feVolLhs->getOpPtrVector().push_back(
1518 // new
1519 // OpCalculateScalarFieldValues(string("FLUXES")+"_residual",resAtGaussPts,MBTET)
1520 // );
1521 feVolLhs->getOpPtrVector().push_back(
1522 new OpValuesAtGaussPts(*this, "VALUES"));
1523 feVolLhs->getOpPtrVector().push_back(
1524 new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
1525 feVolLhs->getOpPtrVector().push_back(
1526 new OpTauDotSigma_HdivHdiv(*this, "FLUXES"));
1527 feVolLhs->getOpPtrVector().push_back(new OpVU_L2L2(*this, "VALUES"));
1528 feVolLhs->getOpPtrVector().push_back(
1529 new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES"));
1530 feVolLhs->getOpPtrVector().push_back(
1531 new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES"));
1532
1533 // Adding finite elements to DM, time solver will ask for it to assemble
1534 // tangent matrices and residuals
1535 boost::shared_ptr<FEMethod> null;
1536 CHKERR DMMoFEMTSSetIFunction(dM, "MIX_BCVALUE", feFaceRhs, null, null);
1537 CHKERR DMMoFEMTSSetIFunction(dM, "MIX", feVolRhs, null, null);
1538 CHKERR DMMoFEMTSSetIJacobian(dM, "MIX", feVolLhs, null, null);
1539
1540 // setting up post-processing
1541
1542
1543 auto get_post_process_ele = [&]() {
1544 auto post_process = boost::make_shared<PostProcEle>(mField);
1545
1546 CHKERR AddHOOps<3, 3, 3>::add(post_process->getOpPtrVector(), {HDIV, L2});
1547
1548 auto values_ptr = boost::make_shared<VectorDouble>();
1549 auto grad_ptr = boost::make_shared<MatrixDouble>();
1550 auto flux_ptr = boost::make_shared<MatrixDouble>();
1551
1552 post_process->getOpPtrVector().push_back(
1553 new OpCalculateScalarFieldValues("VALUES", values_ptr));
1554 post_process->getOpPtrVector().push_back(
1555 new OpCalculateScalarFieldGradient<3>("VALUES", grad_ptr));
1556 post_process->getOpPtrVector().push_back(
1557 new OpCalculateHVecVectorField<3>("FLUXES", flux_ptr));
1558
1559 using OpPPMap = OpPostProcMapInMoab<3, 3>;
1560
1561 post_process->getOpPtrVector().push_back(
1562
1563 new OpPPMap(post_process->getPostProcMesh(),
1564 post_process->getMapGaussPts(),
1565
1566 {{"VALUES", values_ptr}},
1567
1568 {{"GRAD", grad_ptr}, {"FLUXES", flux_ptr}},
1569
1570 {}, {}
1571
1572 ));
1573
1574 return post_process;
1575 };
1576
1577 auto post_process = get_post_process_ele();
1578
1579 post_process->getOpPtrVector().push_back(
1580 new OpValuesAtGaussPts(*this, "VALUES"));
1581 post_process->getOpPtrVector().push_back(
1582 new OpPostProcMaterial(*this, post_process->getPostProcMesh(),
1583 post_process->getMapGaussPts(), "VALUES"));
1584
1585 // Integrate fluxes on boundary
1586 boost::shared_ptr<ForcesAndSourcesCore> flux_integrate;
1587 flux_integrate = boost::shared_ptr<ForcesAndSourcesCore>(
1589 CHKERR AddHOOps<2, 3, 3>::add(flux_integrate->getOpPtrVector(), {HDIV});
1590 flux_integrate->getOpPtrVector().push_back(
1591 new OpIntegrateFluxes(*this, "FLUXES"));
1592 int frequency = 1;
1593 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Monitor post proc", "none");
1594 CHKERR PetscOptionsInt(
1595 "-how_often_output",
1596 "frequency how often results are dumped on hard disk", "", frequency,
1597 &frequency, NULL);
1598 ierr = PetscOptionsEnd();
1599 CHKERRG(ierr);
1600
1601 tsMonitor = boost::shared_ptr<FEMethod>(
1602 new MonitorPostProc(*this, post_process, flux_integrate, frequency));
1603 TsCtx *ts_ctx;
1604 DMMoFEMGetTsCtx(dM, &ts_ctx);
1607 }
1608
1609 Vec D1; ///< Vector with inital head capillary pressure
1610 Vec ghostFlux; ///< Ghost Vector of integrated fluxes
1611
1612 /**
1613 * \brief Create vectors and matrices
1614 * @return Error code
1615 */
1616 MoFEMErrorCode createMatrices() {
1618 CHKERR DMCreateMatrix(dM, &Aij);
1619 CHKERR DMCreateGlobalVector(dM, &D0);
1620 CHKERR VecDuplicate(D0, &D1);
1621 CHKERR VecDuplicate(D0, &D);
1622 CHKERR VecDuplicate(D0, &F);
1623 int ghosts[] = {0};
1624 int nb_locals = mField.get_comm_rank() == 0 ? 1 : 0;
1625 int nb_ghosts = mField.get_comm_rank() > 0 ? 1 : 0;
1626 CHKERR VecCreateGhost(PETSC_COMM_WORLD, nb_locals, 1, nb_ghosts, ghosts,
1627 &ghostFlux);
1629 }
1630
1631 /**
1632 * \brief Delete matrices and vector when no longer needed
1633 * @return error code
1634 */
1635 MoFEMErrorCode destroyMatrices() {
1637 CHKERR MatDestroy(&Aij);
1638 CHKERR VecDestroy(&D);
1639 CHKERR VecDestroy(&D0);
1640 CHKERR VecDestroy(&D1);
1641 CHKERR VecDestroy(&F);
1642 CHKERR VecDestroy(&ghostFlux);
1644 }
1645
1646 /**
1647 * \brief Calculate boundary conditions for fluxes
1648 * @return Error code
1649 */
1650 MoFEMErrorCode calculateEssentialBc() {
1652 // clear vectors
1653 CHKERR VecZeroEntries(D0);
1654 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
1655 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
1656 // clear essential bc indices, it could have dofs from other mesh refinement
1657 bcIndices.clear();
1658 // set operator to calculate essential boundary conditions
1659 CHKERR DMoFEMLoopFiniteElements(dM, "MIX_BCFLUX", feFaceBc);
1660 CHKERR VecAssemblyBegin(D0);
1661 CHKERR VecAssemblyEnd(D0);
1662 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
1663 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
1664 double norm2D0;
1665 CHKERR VecNorm(D0, NORM_2, &norm2D0);
1666 // CHKERR VecView(D0,PETSC_VIEWER_STDOUT_WORLD);
1667 PetscPrintf(PETSC_COMM_WORLD, "norm2D0 = %6.4e\n", norm2D0);
1669 }
1670
1671 /**
1672 * \brief Calculate inital pressure head distribution
1673 * @return Error code
1674 */
1675 MoFEMErrorCode calculateInitialPc() {
1677 // clear vectors
1678 CHKERR VecZeroEntries(D1);
1679 CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_FORWARD);
1680 CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_FORWARD);
1681 // Calculate initial pressure head on each element
1682 CHKERR DMoFEMLoopFiniteElements(dM, "MIX", feVolInitialPc);
1683 // Assemble vector
1684 CHKERR VecAssemblyBegin(D1);
1685 CHKERR VecAssemblyEnd(D1);
1686 CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_FORWARD);
1687 CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_FORWARD);
1688 // Calculate norm
1689 double norm2D1;
1690 CHKERR VecNorm(D1, NORM_2, &norm2D1);
1691 // CHKERR VecView(D0,PETSC_VIEWER_STDOUT_WORLD);
1692 PetscPrintf(PETSC_COMM_WORLD, "norm2D1 = %6.4e\n", norm2D1);
1694 }
1695
1696 /**
1697 * \brief solve problem
1698 * @return error code
1699 */
1700 MoFEMErrorCode solveProblem(bool set_initial_pc = true) {
1702 if (set_initial_pc) {
1703 // Set initial head
1704 CHKERR DMoFEMMeshToLocalVector(dM, D1, INSERT_VALUES, SCATTER_REVERSE);
1705 }
1706
1707 // Initiate vector from data on the mesh
1708 CHKERR DMoFEMMeshToLocalVector(dM, D, INSERT_VALUES, SCATTER_FORWARD);
1709
1710 // Create time solver
1711 TS ts;
1712 CHKERR TSCreate(PETSC_COMM_WORLD, &ts);
1713 // Use backward Euler method
1714 CHKERR TSSetType(ts, TSBEULER);
1715 // Set final time
1716 double ftime = 1;
1717 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
1718 // Setup solver from commabd line
1719 CHKERR TSSetFromOptions(ts);
1720 // Set DM to TS
1721 CHKERR TSSetDM(ts, dM);
1722#if PETSC_VERSION_GE(3, 7, 0)
1723 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
1724#endif
1725 // Set-up monitor
1726 TsCtx *ts_ctx;
1727 DMMoFEMGetTsCtx(dM, &ts_ctx);
1728 CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx, PETSC_NULL);
1729
1730 // This add SNES monitor, to show error by fields. It is dirty trick
1731 // to add monitor, so code is hidden from doxygen
1732 CHKERR TSSetSolution(ts, D);
1733 CHKERR TSSetUp(ts);
1734 SNES snes;
1735 CHKERR TSGetSNES(ts, &snes);
1736
1737#if PETSC_VERSION_GE(3, 7, 0)
1738 {
1739 PetscViewerAndFormat *vf;
1740 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
1741 PETSC_VIEWER_DEFAULT, &vf);
1742 CHKERR SNESMonitorSet(
1743 snes,
1744 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
1745 void *))SNESMonitorFields,
1746 vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
1747 }
1748#else
1749 {
1750 CHKERR SNESMonitorSet(snes,
1751 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
1752 void *))SNESMonitorFields,
1753 0, 0);
1754 }
1755#endif
1756
1757 CHKERR TSSolve(ts, D);
1758
1759 // Get statistic form TS and print it
1760 CHKERR TSGetTime(ts, &ftime);
1761 PetscInt steps, snesfails, rejects, nonlinits, linits;
1762 CHKERR TSGetTimeStepNumber(ts, &steps);
1763 CHKERR TSGetSNESFailures(ts, &snesfails);
1764 CHKERR TSGetStepRejections(ts, &rejects);
1765 CHKERR TSGetSNESIterations(ts, &nonlinits);
1766 CHKERR TSGetKSPIterations(ts, &linits);
1767 PetscPrintf(PETSC_COMM_WORLD,
1768 "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits "
1769 "%D, linits %D\n",
1770 steps, rejects, snesfails, ftime, nonlinits, linits);
1771
1773 }
1774};
1775
1776} // namespace MixTransport
1777
1778#endif // __UNSATURATD_FLOW_HPP__
static Index< 'K', 3 > K
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
static PetscErrorCode ierr
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
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 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 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_row(const std::string &fe_name, const std::string name_row)=0
set field row 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
double D
const double v
phase velocity of light in medium (cm/ns)
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1876
PostProcBrokenMeshInMoab< VolumeElementForcesAndSourcesCore > PostProcEle
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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 pressure 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.
BasicMethodsSequence & getPostProcessMonitor()
Get the postProcess to do Monitor object.
Definition: TsCtx.hpp:154
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Force scale operator for reading two columns.