v0.14.0
Searching...
No Matches
Gels.hpp
Go to the documentation of this file.
1/** \file Gels.hpp
2 \brief Implementation of Gel finite element
3 \ingroup gel
4*/
5
6/* This file is part of MoFEM.
7 * MoFEM is free software: you can redistribute it and/or modify it under
9 * Free Software Foundation, either version 3 of the License, or (at your
10 * option) any later version.
11 *
12 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
13 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15 * License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public
19
20#ifndef __GEL_HPP__
21#define __GEL_HPP__
22
24 #error "MoFEM need to be compiled with ADOL-C"
25#endif
26
27namespace GelModule {
28
29/** \brief Implementation of Gel constitutive model
30\ingroup gel
31
32Implementation follows constitutive equation from:
33
34VISCOELASTICITY AND POROELASTICITY IN ELASTOMERIC GELS
35Acta Mechanica Solida Sinica, Vol. 25, No. 5,
36Yuhang Hu Zhigang Suo
37
38Note1: Following implementation is tailored for large strain analysis,
39however in current version the engineering strain is used, i.e. assuming small deformation.
40Moreover calculation of solvent flux and other physical quantities is with
41assumption that those values are nonlinear.
42
43Note2: Basic implementation is with linear constitutive equations, however
44adding nonlinearities to model should be simple, since all tangent matrices
46
47System of equations of equations in the form is assembled:
48\f[
49\left[
50\begin{array}{ccc}
51\mathbf{K}_{xx} & \mathbf{K}_{x\mu} & \mathbf{K}_{x\hat{\varepsilon}} \\
52\mathbf{K}_{\mu x} & \mathbf{K}_{\mu \mu} & \mathbf{K}_{\mu \hat{\varepsilon}} \\
53\mathbf{K}_{\hat{\varepsilon} x} & \mathbf{K}_{\hat{\varepsilon} \mu} & \mathbf{K}_{\hat{\varepsilon} \hat{\varepsilon}}
54\end{array}
55\right]
56\left[
57\begin{array}{c}
58\mathbf{q}_{x}\\
59\mathbf{q}_{\mu}\\
60\mathbf{q}_{\hat{\varepsilon}}
61\end{array}
62\right]
63=
64\left[
65\begin{array}{l}
66\mathbf{f}^\textrm{internal}_{x}\\
67\mathbf{f}^V_{\mu}+\mathbf{f}^J_{\mu}\\
68\mathbf{r}_{\hat{\varepsilon}}
69\end{array}
70\right]
71\f]
72
73*/
74struct Gel {
75
81 };
82
84
85 /** \brief Gel material parameters
86 \ingroup gel
87 */
89
90 // FIXME: Need to have hood description for each parameter and units
91 double vAlpha; ///< Poisson ratio spring alpha
92 double gAlpha; ///< Shear modulus spring alpha
93 double vBeta; ///< Poisson ratio spring beta
94 double gBeta; ///< Shear modulus spring beta
95 double vBetaHat; ///< Poisson ratio dashpot
96 double gBetaHat; ///< Shear modulus dashpot
97 double pErmeability; ///< Permeability of Gel
98 double vIscosity; ///< Viscosity of Gel
99 double oMega; ///< Volume per solvent molecule
100 double mU0; ///< equilibrated chemical potential/load
101
102 Range tEts; ///< Test in the block
103
104 };
105 map<int,BlockMaterialData> blockMaterialData;
106
107 /** \brief Constitutive model functions
108 \ingroup gel
109
110 \image html gel_spring_daspot_model.png "Gel model" width=700px
111
112 */
113 template<typename TYPE>
115
116 int iD;
117 map<int,BlockMaterialData> &dAta;
118
119 ConstitutiveEquation(map<int,BlockMaterialData> &data):
120 dAta(data) {
121 }
123
124 // Input
125
126 ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > F; ///< Gradient of deformation
127 ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > FDot; ///< Rate of gradient of deformation
128 ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > strainHat; ///< Internal variable, strain in dashpot beta
129 ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > strainHatDot; ///< Internal variable, strain in dashpot beta
130 TYPE mU; ///< Solvent concentration
132
133 // Internal use
134
135 ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > C; ///< Cauchy deformation
137 ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > strainTotal; ///< Total strain applied at integration point
138
139 ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > stressAlpha; ///< Stress generated by spring alpha
140 ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > stressBeta; ///< Stress generated by spring beta
141 ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > strainHatFlux; ///< Rate of dashpot (beta) strain
142 ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > stressBetaHat; ///< Stress as result of volume change due to solvent concentration
143
148
149 // Output
150
151 ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > stressTotal; ///< Total stress
152 ublas::vector<TYPE,ublas::bounded_array<TYPE,9> > solventFlux; ///< Solvent flux
153 TYPE solventConcentrationDot; ///< Volume rate change
154 ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > residualStrainHat; ///< Residual for calculation epsilon hat
155
156 virtual PetscErrorCode calculateCauchyDefromationTensor() {
157 PetscFunctionBegin;
158 C.resize(3,3);
159 noalias(C) = prod(trans(F),F);
160 PetscFunctionReturn(0);
161 }
162
163 /** \brief Calculate total strain
164
165 */
166 virtual PetscErrorCode calculateStrainTotal() {
167 PetscFunctionBegin;
170 for(int ii = 0;ii<3;ii++) {
172 }
173 strainTotal.resize(3,3,false);
175 strainTotal *= 0.5;
176 PetscFunctionReturn(0);
177 }
178
179 /** \berief Calculate trace of rate of gradient of deformation
180 */
181 virtual PetscErrorCode calculateTraceStrainTotalDot() {
182 PetscFunctionBegin;
184 for(int ii = 0;ii<3;ii++) {
185 traceStrainTotalDot += FDot(ii,ii);
186 }
187 PetscFunctionReturn(0);
188 }
189
190 /** \brief Calculate stress in spring alpha
191
192 \f[
193 \sigma^\alpha_{ij} = 2G^\alpha\left(\varepsilon_{ij} + \frac{v^\alpha}{1-2v^\alpha}\varepsilon_{kk}\delta_{ij}\right)
194 \f]
195
196 Note: In general implementation for large strain this function should calculate
197 Piola-Kirchhoff Stress I.
198
199 */
200 virtual PetscErrorCode calculateStressAlpha() {
201 PetscFunctionBegin;
203 for(int ii = 0;ii<3;ii++) {
205 }
206 stressAlpha.resize(3,3,false);
207 double a = 2.0*dAta[iD].gAlpha;
208 double b = a*(dAta[iD].vAlpha/(1.0-2.0*dAta[iD].vAlpha));
209 noalias(stressAlpha) = a*strainTotal;
210 for(int ii = 0; ii<3; ii++){
211 stressAlpha(ii,ii) += b*traceStrainTotal;
212 }
213 PetscFunctionReturn(0);
214 }
215
216 /** \brief Calculate stress in spring beta
217
218 \f[
219 \sigma^\beta_{ij} = 2G^\beta\left[
220 (\varepsilon_{ij}-\hat{\varepsilon}_{ij})
221 + \frac{v^\beta}{1-2v^\beta}(\varepsilon_{kk}-\hat{\varepsilon}_{kk})\delta_{ij}
222 \right]
223 \f]
224
225 */
226 virtual PetscErrorCode calculateStressBeta() {
227 PetscFunctionBegin;
229 traceStrainHat = 0;
230 for(int ii = 0;ii<3;ii++) {
231 traceStrainHat += strainHat(ii,ii);
233 }
234 stressBeta.resize(3,3,false);
235 double a = 2.0*dAta[iD].gBeta;
236 double b = a*(dAta[iD].vBeta/(1.0-2.0*dAta[iD].vBeta));
237 noalias(stressBeta) = a*(strainTotal-strainHat);
238 for(int ii = 0;ii<3;ii++) {
240 }
241 PetscFunctionReturn(0);
242 }
243
244 /** \brief Calculate rate of strain hat
245
246 \f[
247 \frac{\partial \hat{\varepsilon}_{ij}}{\partial t} =
248 \frac{1}{2\hat{G}^\beta}
249 \left(
250 \sigma_{ij}^\beta - \frac{\hat{v}^\beta}{1+\hat{v}^\beta}\sigma_{kk}^\beta\delta_{ij}
251 \right)
252 \f]
253
254 */
255 virtual PetscErrorCode calculateStrainHatFlux() {
256 PetscFunctionBegin;
257 traceStressBeta = 0;
258 for(int ii = 0;ii<3;ii++) {
260 }
261 strainHatFlux.resize(3,3,false);
262 double a = (1.0/(2.0*dAta[iD].gBetaHat));
263 noalias(strainHatFlux) = a*stressBeta;
264 double b = a*(dAta[iD].vBetaHat/(1.0+dAta[iD].vBetaHat));
265 for(int ii = 0;ii<3;ii++) {
266 strainHatFlux(ii,ii) -= b*traceStressBeta;
267 }
268 PetscFunctionReturn(0);
269 }
270
271 /** \brief Calculate stress due to concentration of solvent molecules
272
273 \f[
274 \hat{\sigma}_{ij}^\beta = \frac{\Delta\mu}{\Omega}\delta_{ij}
275 \f]
276 note that \f$\Delta\mu = \mu-\mu_0\f$. In current version of code
277 value of \f$\Delta\mu\f$ is approximated.
278
279 */
280 virtual PetscErrorCode calculateStressBetaHat() {
281 PetscFunctionBegin;
282 stressBetaHat.resize(3,3,false);
283 stressBetaHat.clear();
284 for(int ii=0; ii<3; ii++){
285 stressBetaHat(ii,ii) = -(mU-dAta[iD].mU0)/dAta[iD].oMega;
286 }
287 PetscFunctionReturn(0);
288 }
289
290 // Functions calculating output variables
291
292 virtual PetscErrorCode calculateStressTotal() {
293 PetscFunctionBegin;
294 stressTotal.resize(3,3,false);
295 noalias(stressTotal) = stressAlpha;
296 noalias(stressTotal) += stressBeta;
297 noalias(stressTotal) += stressBetaHat;
298 PetscFunctionReturn(0);
299 }
300
301 virtual PetscErrorCode calculateResidualStrainHat() {
302 PetscFunctionBegin;
303 residualStrainHat.resize(3,3,false);
305 PetscFunctionReturn(0);
306 }
307
308 /** \brief Calculate flux
309
310 \f[
311 J_k =
312 -\left(
313 \frac{\kappa}{\eta\Omega^2}
314 \right)
315 \mu_{,i}
316 \f]
317
318 */
319 virtual PetscErrorCode calculateSolventFlux() {
320 PetscFunctionBegin;
321 double a = dAta[iD].pErmeability/(dAta[iD].vIscosity*dAta[iD].oMega/**dAta[iD].oMega*/);
323 PetscFunctionReturn(0);
324 }
325
326 /** \brief Calculate solvent concentration rate
327
328 FIXME: For simplicity as first approximation set volume change
329 as trace of gradient total strain
330
331 */
332 virtual PetscErrorCode calculateSolventConcentrationDot() {
333 PetscFunctionBegin;
334 solventConcentrationDot = traceStrainTotalDot/*/dAta[iD].oMega*/;
335 PetscFunctionReturn(0);
336 }
337
338 };
339
341
342 /** \brief Common data for gel model
343 \ingroup gel
344 */
345 struct CommonData {
346
351 string muName;
352 string muNameDot;
353
354 map<string,vector<VectorDouble> > dataAtGaussPts;
356
357 vector<MatrixDouble3by3> stressTotal;
358 vector<VectorDouble3> solventFlux;
360 vector<VectorDouble9> residualStrainHat;
361
362 vector<double*> jacRowPtr;
363 vector<MatrixDouble> jacStressTotal;
364 vector<MatrixDouble> jacSolventFlux;
365 vector<VectorDouble> jacSolventDot;
366 vector<MatrixDouble> jacStrainHat;
367
370
372 recordOn(true) {
373 }
374
375 };
377
378 /// \brief definition of volume element
380
382 int addToRule; ///< Takes into account HO geometry
383
384 GelFE(MoFEM::Interface &m_field,CommonData &common_data):
386 commonData(common_data),
388 }
389
390 int getRule(int order) {
392 }
393
394 PetscErrorCode preProcess() {
395
396 PetscFunctionBegin;
398
400
401 ierr = mField.set_other_local_ghost_vector(
403 ); CHKERRQ(ierr);
404 ierr = mField.set_other_local_ghost_vector(
405 problemPtr,commonData.muName,commonData.muNameDot,COL,ts_u_t,INSERT_VALUES,SCATTER_REVERSE
406 ); CHKERRQ(ierr);
407 ierr = mField.set_other_local_ghost_vector(
409 ); CHKERRQ(ierr);
410
411 }
412
413 PetscFunctionReturn(0);
414 }
415
416 PetscErrorCode postProcess() {
417
418 PetscFunctionBegin;
419
421
423 ierr = VecAssemblyBegin(ts_F); CHKERRQ(ierr);
424 ierr = VecAssemblyEnd(ts_F); CHKERRQ(ierr);
425 }
427 ierr = MatAssemblyBegin(ts_B,MAT_FLUSH_ASSEMBLY); CHKERRQ(ierr);
428 ierr = MatAssemblyEnd(ts_B,MAT_FLUSH_ASSEMBLY); CHKERRQ(ierr);
429 }
430
431 PetscFunctionReturn(0);
432 }
433
434 };
435
437
439 mFiled(m_field),
440 feRhs(m_field,commonData),
441 feLhs(m_field,commonData) {
442 }
443
445
450
452 const string field_name,
453 CommonData &common_data,
454 bool calc_val,
456 EntityType zero_at_type = MBVERTEX
457 ):
458 MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
460 ),
461 commonData(common_data),
462 calcVal(calc_val),
464 zeroAtType(zero_at_type) {
465 }
466
467 /** \brief Operator field value
468 *
469 */
470 PetscErrorCode doWork(
471 int side,EntityType type,DataForcesAndSourcesCore::EntData &data
472 ){
473 PetscFunctionBegin;
474
475 try {
476
477 int nb_dofs = data.getFieldData().size();
478 if(nb_dofs == 0) {
479 PetscFunctionReturn(0);
480 }
481 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
482 int nb_gauss_pts = data.getN().size1();
483
484 // Initialize
485 if(calcVal) {
486 commonData.dataAtGaussPts[rowFieldName].resize(nb_gauss_pts);
487 for(int gg = 0;gg<nb_gauss_pts;gg++) {
488 commonData.dataAtGaussPts[rowFieldName][gg].resize(rank,false);
489 }
490 }
493 for(int gg = 0;gg<nb_gauss_pts;gg++) {
495 }
496 }
497
498 // Zero values
499 if(type == zeroAtType) {
500 for(int gg = 0;gg<nb_gauss_pts;gg++) {
501 if(calcVal) {
503 }
506 }
507 }
508 }
509
510 VectorDouble &values = data.getFieldData();
511
512 // Calculate values at integration points
513 for(int gg = 0;gg<nb_gauss_pts;gg++) {
514 VectorDouble N = data.getN(gg,nb_dofs/rank);
515 MatrixDouble diffN = data.getDiffN(gg,nb_dofs/rank);
516 for(int dd = 0;dd<nb_dofs/rank;dd++) {
517 for(int rr1 = 0;rr1<rank;rr1++) {
518 if(calcVal) {
519 commonData.dataAtGaussPts[rowFieldName][gg][rr1] += N[dd]*values[rank*dd+rr1];
520 }
522 for(int rr2 = 0;rr2<3;rr2++) {
524 }
525 }
526 }
527 }
528 }
529
530 } catch (const std::exception& ex) {
531 ostringstream ss;
532 ss << "throw in method: " << ex.what() << endl;
533 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
534 }
535
536 PetscFunctionReturn(0);
537 }
538
539 };
540
542
543 vector<int> tagS;
546
549 bool &recordOn;
550 map<int,int> &nbActiveVariables;
551 map<int,int> &nbActiveResults;
552
554 const string field_name,
555 vector<int> tags,
557 CommonData &common_data,
558 bool calculate_residual,
559 bool calculate_jacobian
560 ):
561 MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
563 ),
564 tagS(tags),
565 cE(ce),
566 commonData(common_data),
567 calculateResidualBool(calculate_residual),
568 calculateJacobianBool(calculate_jacobian),
569 recordOn(common_data.recordOn),
571 nbActiveResults(common_data.nbActiveResults) {
572 }
573
575 VectorDouble activeVariables;
576
577 PetscErrorCode recordStressTotal() {
578
579 PetscFunctionBegin;
580
581 try {
582 if(tagS[STRESSTOTAL]<0) {
583 PetscFunctionReturn(0);
584 }
585 cE->F.resize(3,3,false);
586 cE->strainHat.resize(3,3,false);
588 VectorDouble &strain_hat = (commonData.dataAtGaussPts[commonData.strainHatName])[0];
589 VectorDouble mu = (commonData.dataAtGaussPts[commonData.muName])[0];
590 trace_on(tagS[STRESSTOTAL]);
591 {
592 // Activate gradient of defamation
594 for(int dd1 = 0;dd1<3;dd1++) {
595 for(int dd2 = 0;dd2<3;dd2++) {
596 cE->F(dd1,dd2) <<= F(dd1,dd2);
598 }
599 }
600 // Using vector notation to store hatStrain
601 // xx,yy,zz,2xy,2yz,2zy
602 cE->strainHat(0,0) <<= strain_hat[0];
603 cE->strainHat(1,1) <<= strain_hat[1];
604 cE->strainHat(2,2) <<= strain_hat[2];
605 cE->strainHat(0,1) <<= strain_hat[3];
606 cE->strainHat(1,2) <<= strain_hat[4];
607 cE->strainHat(0,2) <<= strain_hat[5];
609 cE->strainHat(1,0) = cE->strainHat(0,1);
610 cE->strainHat(2,1) = cE->strainHat(1,2);
611 cE->strainHat(2,0) = cE->strainHat(0,2);
612 cE->mU <<= mu[0];
614 // Do calculations
615 ierr = cE->calculateStrainTotal(); CHKERRQ(ierr);
616 ierr = cE->calculateStressAlpha(); CHKERRQ(ierr);
617 ierr = cE->calculateStressBeta(); CHKERRQ(ierr);
618 ierr = cE->calculateStressBetaHat(); CHKERRQ(ierr);
619 // Finally calculate result
620 ierr = cE->calculateStressTotal(); CHKERRQ(ierr);
621 // Results
624 commonData.stressTotal[0].resize(3,3,false);
625 for(int d1 = 0;d1<3;d1++) {
626 for(int d2 = 0;d2<3;d2++) {
627 cE->stressTotal(d1,d2) >>= (commonData.stressTotal[0])(d1,d2);
629 }
630 }
631 }
632 trace_off();
633 } catch (const std::exception& ex) {
634 ostringstream ss;
635 ss << "throw in method: " << ex.what() << endl;
636 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
637 }
638 PetscFunctionReturn(0);
639 }
640
641 PetscErrorCode recordSolventFlux() {
642 PetscFunctionBegin;
643 try {
644 if(tagS[SOLVENTFLUX]<0) {
645 PetscFunctionReturn(0);
646 }
647
650 trace_on(tagS[SOLVENTFLUX]);
651 {
652 // Activate rate of gradient of defamation
654 for(int ii = 0;ii<3;ii++) {
657 }
658 ierr = cE->calculateSolventFlux(); CHKERRQ(ierr);
661 commonData.solventFlux[0].resize(3,false);
662 for(int d1 = 0;d1<3;d1++) {
663 cE->solventFlux[d1] >>= commonData.solventFlux[0][d1];
665 }
666 }
667 trace_off();
668 } catch (const std::exception& ex) {
669 ostringstream ss;
670 ss << "throw in method: " << ex.what() << endl;
671 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
672 }
673
674 PetscFunctionReturn(0);
675 }
676
678 PetscFunctionBegin;
679
680 if(tagS[SOLVENTRATE]<0) {
681 PetscFunctionReturn(0);
682 }
683
684
685
688
689 trace_on(tagS[SOLVENTRATE]);
690 {
691 // Activate rate of gradient of defamation
693 cE->F.resize(3,3,false);
694 for(int dd1 = 0;dd1<3;dd1++) {
695 for(int dd2 = 0;dd2<3;dd2++) {
696 cE->F(dd1,dd2) <<= F(dd1,dd2);
698 }
699 }
700 cE->FDot.resize(3,3,false);
701 for(int dd1 = 0;dd1<3;dd1++) {
702 for(int dd2 = 0;dd2<3;dd2++) {
703 cE->FDot(dd1,dd2) <<= F_dot(dd1,dd2);
705 }
706 }
707 ierr = cE->calculateTraceStrainTotalDot(); CHKERRQ(ierr);
708 ierr = cE->calculateSolventConcentrationDot(); CHKERRQ(ierr);
711 cE->solventConcentrationDot >>= commonData.solventConcentrationDot[0];
713 }
714 trace_off();
715
716 PetscFunctionReturn(0);
717 }
718
719 PetscErrorCode recordResidualStrainHat() {
720 PetscFunctionBegin;
721
722 if(tagS[RESIDUALSTRAINHAT]<0) {
723 PetscFunctionReturn(0);
724 }
725
726
727
728 cE->F.resize(3,3,false);
729 cE->strainHat.resize(3,3,false);
730 cE->strainHatDot.resize(3,3,false);
732 VectorDouble &strain_hat = (commonData.dataAtGaussPts[commonData.strainHatName])[0];
733 VectorDouble &strain_hat_dot = (commonData.dataAtGaussPts[commonData.strainHatNameDot])[0];
734 trace_on(tagS[RESIDUALSTRAINHAT]);
735 {
736 // Activate gradient of defamation
738 for(int dd1 = 0;dd1<3;dd1++) {
739 for(int dd2 = 0;dd2<3;dd2++) {
740 cE->F(dd1,dd2) <<= F(dd1,dd2);
742 }
743 }
744 // Using vector notation to store hatStrain
745 // xx,yy,zz,2xy,2yz,2zy
746 cE->strainHat(0,0) <<= strain_hat[0];
747 cE->strainHat(1,1) <<= strain_hat[1];
748 cE->strainHat(2,2) <<= strain_hat[2];
749 cE->strainHat(0,1) <<= strain_hat[3];
750 cE->strainHat(1,2) <<= strain_hat[4];
751 cE->strainHat(0,2) <<= strain_hat[5];
753 cE->strainHat(1,0) = cE->strainHat(0,1);
754 cE->strainHat(2,1) = cE->strainHat(1,2);
755 cE->strainHat(2,0) = cE->strainHat(0,2);
756 // Activate strain hat dot
757 cE->strainHatDot(0,0) <<= strain_hat_dot[0];
758 cE->strainHatDot(1,1) <<= strain_hat_dot[1];
759 cE->strainHatDot(2,2) <<= strain_hat_dot[2];
760 cE->strainHatDot(0,1) <<= strain_hat_dot[3];
761 cE->strainHatDot(1,2) <<= strain_hat_dot[4];
762 cE->strainHatDot(0,2) <<= strain_hat_dot[5];
764 cE->strainHatDot(1,0) = cE->strainHatDot(0,1);
765 cE->strainHatDot(2,1) = cE->strainHatDot(1,2);
766 cE->strainHatDot(2,0) = cE->strainHatDot(0,2);
767 ierr = cE->calculateStrainTotal(); CHKERRQ(ierr);
768 ierr = cE->calculateStressBeta(); CHKERRQ(ierr);
769 ierr = cE->calculateStrainHatFlux(); CHKERRQ(ierr);
770 ierr = cE->calculateResidualStrainHat(); CHKERRQ(ierr);
773 commonData.residualStrainHat[0].resize(6,false);
774 cE->residualStrainHat(0,0) >>= commonData.residualStrainHat[0][nbActiveResults[tagS[RESIDUALSTRAINHAT]]++];
775 cE->residualStrainHat(1,1) >>= commonData.residualStrainHat[0][nbActiveResults[tagS[RESIDUALSTRAINHAT]]++];
776 cE->residualStrainHat(2,2) >>= commonData.residualStrainHat[0][nbActiveResults[tagS[RESIDUALSTRAINHAT]]++];
777 cE->residualStrainHat(0,1) >>= commonData.residualStrainHat[0][nbActiveResults[tagS[RESIDUALSTRAINHAT]]++];
778 cE->residualStrainHat(1,2) >>= commonData.residualStrainHat[0][nbActiveResults[tagS[RESIDUALSTRAINHAT]]++];
779 cE->residualStrainHat(0,2) >>= commonData.residualStrainHat[0][nbActiveResults[tagS[RESIDUALSTRAINHAT]]++];
780 }
781 trace_off();
782
783 PetscFunctionReturn(0);
784 }
785
786 PetscErrorCode calculateFunction(TagEvaluate te,double *ptr) {
787 PetscFunctionBegin;
788
789 int r;
790 //play recorder for values
791 r = ::function(
792 tagS[te],
795 &activeVariables[0],
796 ptr
797 );
798 if(r<3) { // function is locally analytic
799 SETERRQ1(
800 PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,
801 "ADOL-C function evaluation with error r = %d",r
802 );
803 }
804
805 PetscFunctionReturn(0);
806 }
807
808 PetscErrorCode calculateJacobian(TagEvaluate te) {
809 PetscFunctionBegin;
810 try {
811 int r;
812 r = jacobian(
813 tagS[te],
816 &activeVariables[0],
818 );
819 if(r<3) {
820 SETERRQ(
821 PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,
822 "ADOL-C function evaluation with error"
823 );
824 }
825 } catch (const std::exception& ex) {
826 ostringstream ss;
827 ss << "throw in method: " << ex.what() << endl;
828 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
829 }
830 PetscFunctionReturn(0);
831 }
832
834 PetscFunctionBegin;
835
836 try {
837
838 if(tagS[STRESSTOTAL]<0) {
839 PetscFunctionReturn(0);
840 }
841
842
843
845
846 for(int gg = 0;gg<nbGaussPts;gg++) {
847
849 VectorDouble &strain_hat = (commonData.dataAtGaussPts[commonData.strainHatName])[gg];
850 double mu = (commonData.dataAtGaussPts[commonData.muName])[gg][0];
851 int nb_active_variables = 0;
852 // Activate gradient of defamation
853 for(int dd1 = 0;dd1<3;dd1++) {
854 for(int dd2 = 0;dd2<3;dd2++) {
855 activeVariables[nb_active_variables++] = F(dd1,dd2);
856 }
857 }
858 // Using vector notation to store hatStrain
859 // xx,yy,zz,2xy,2yz,2zy
860 for(int ii = 0;ii<6;ii++) {
861 activeVariables[nb_active_variables++] = strain_hat[ii];
862 }
864 activeVariables[nb_active_variables++] = mu;
865
866 if(nb_active_variables!=nbActiveVariables[tagS[STRESSTOTAL]]) {
867 SETERRQ(
868 PETSC_COMM_SELF,MOFEM_IMPOSIBLE_CASE,
869 "Number of active variables does not much"
870 );
871 }
872
874 if(gg == 0) {
876 }
877 commonData.stressTotal[gg].resize(3,3,false);
879 }
880
882 if(gg == 0) {
885 }
886 commonData.jacStressTotal[gg].resize(
889 false
890 );
891 for(int dd = 0;dd<nbActiveResults[tagS[STRESSTOTAL]];dd++) {
893 }
895 }
896
897 }
898
899 } catch (const std::exception& ex) {
900 ostringstream ss;
901 ss << "throw in method: " << ex.what() << endl;
902 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
903 }
904
905 PetscFunctionReturn(0);
906 }
907
909 PetscFunctionBegin;
910
911 if(tagS[SOLVENTFLUX]<0) {
912 PetscFunctionReturn(0);
913 }
914
915
916
918
919 for(int gg = 0;gg<nbGaussPts;gg++) {
921 int nb_active_variables = 0;
922 // Activate rate of gradient of defamation
923 for(int ii = 0;ii<3;ii++) {
925 }
926 if(nb_active_variables!=nbActiveVariables[tagS[SOLVENTFLUX]]) {
927 SETERRQ(
928 PETSC_COMM_SELF,MOFEM_IMPOSIBLE_CASE,"Number of active variables does not much"
929 );
930 }
932 commonData.solventFlux[gg].resize(3,false);
935 &(commonData.solventFlux[gg][0])
936 ); CHKERRQ(ierr);
937 }
939 if(gg == 0) {
942 }
943 commonData.jacSolventFlux[gg].resize(
946 false
947 );
948 for(int dd = 0;dd<nbActiveResults[tagS[SOLVENTFLUX]];dd++) {
950 }
952 }
953 }
954
955 PetscFunctionReturn(0);
956 }
957
959 PetscFunctionBegin;
960
961 if(tagS[SOLVENTRATE]<0) {
962 PetscFunctionReturn(0);
963 }
964
965
966
968
969 for(int gg = 0;gg<nbGaussPts;gg++) {
972 int nb_active_variables = 0;
973 // Activate gradient of defamation
974 for(int dd1 = 0;dd1<3;dd1++) {
975 for(int dd2 = 0;dd2<3;dd2++) {
976 activeVariables[nb_active_variables++] = F(dd1,dd2);
977 }
978 }
979 for(int dd1 = 0;dd1<3;dd1++) {
980 for(int dd2 = 0;dd2<3;dd2++) {
981 activeVariables[nb_active_variables++] = F_dot(dd1,dd2);
982 }
983 }
984 if(nb_active_variables!=nbActiveVariables[tagS[SOLVENTRATE]]) {
985 SETERRQ(
986 PETSC_COMM_SELF,MOFEM_IMPOSIBLE_CASE,"Number of active variables does not much"
987 );
988 }
991 }
993 if(gg == 0) {
996 }
1000 }
1001 }
1002
1003 PetscFunctionReturn(0);
1004 }
1005
1007 PetscFunctionBegin;
1008
1009 try {
1010 if(tagS[RESIDUALSTRAINHAT]<0) {
1011 PetscFunctionReturn(0);
1012 }
1013
1015 for(int gg = 0;gg<nbGaussPts;gg++) {
1017 VectorDouble &strain_hat = (commonData.dataAtGaussPts[commonData.strainHatName])[gg];
1018 VectorDouble &strain_hat_dot = (commonData.dataAtGaussPts[commonData.strainHatNameDot])[gg];
1019 int nb_active_variables = 0;
1020 // Activate gradient of defamation
1021 for(int dd1 = 0;dd1<3;dd1++) {
1022 for(int dd2 = 0;dd2<3;dd2++) {
1023 activeVariables[nb_active_variables++] = F(dd1,dd2);
1024 }
1025 }
1026 // Using vector notation to store hatStrain
1027 // xx,yy,zz,2xy,2yz,2zy
1028 for(int ii = 0;ii<6;ii++) {
1029 activeVariables[nb_active_variables++] = strain_hat[ii];
1030 }
1031 // Using vector notation to store hatStrainDot
1032 for(int ii = 0;ii<6;ii++) {
1033 activeVariables[nb_active_variables++] = strain_hat_dot[ii];
1034 }
1035 if(nb_active_variables!=nbActiveVariables[tagS[RESIDUALSTRAINHAT]]) {
1036 SETERRQ(
1037 PETSC_COMM_SELF,MOFEM_IMPOSIBLE_CASE,"Number of active variables does not much"
1038 );
1039 }
1041 if(gg == 0) {
1043 }
1044 commonData.residualStrainHat[gg].resize(6,false);
1047 ); CHKERRQ(ierr);
1048 }
1050 if(gg == 0) {
1053 }
1054 commonData.jacStrainHat[gg].resize(
1057 false
1058 );
1059 for(int dd = 0;dd<nbActiveResults[tagS[RESIDUALSTRAINHAT]];dd++) {
1061 }
1063 }
1064 }
1065
1066
1067 } catch (const std::exception& ex) {
1068 ostringstream ss;
1069 ss << "throw in method: " << ex.what() << endl;
1070 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1071 }
1072
1073 PetscFunctionReturn(0);
1074 }
1075
1076 PetscErrorCode doWork(
1077 int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
1078 ) {
1079 PetscFunctionBegin;
1080
1081
1082
1083 if(row_type != MBVERTEX) PetscFunctionReturn(0);
1084 nbGaussPts = row_data.getN().size1();
1085
1086 Tag th_block_id;
1088 rval = getVolumeFE()->
1089 mField.get_moab().tag_get_handle("BLOCK_ID",th_block_id); CHKERRQ_MOAB(rval);
1090 rval = getVolumeFE()->
1091 mField.get_moab().tag_get_data(th_block_id,&ent,1,&(cE->iD)); CHKERRQ_MOAB(rval);
1092
1093 try {
1094
1095 if(recordOn) {
1096 ierr = recordStressTotal(); CHKERRQ(ierr);
1097 ierr = recordSolventFlux(); CHKERRQ(ierr);
1099 ierr = recordResidualStrainHat(); CHKERRQ(ierr);
1100 }
1101
1106
1107 } catch (const std::exception& ex) {
1108 ostringstream ss;
1109 ss << "throw in method: " << ex.what() << endl;
1110 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1111 }
1112
1113 PetscFunctionReturn(0);
1114 }
1115
1116 };
1117
1120 MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(field_name,UserDataOperator::OPROW) {
1121 }
1122
1123 VectorDouble nF;
1124 PetscErrorCode aSemble(
1125 int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
1126 ) {
1127 PetscFunctionBegin;
1128 try {
1129 int nb_dofs = row_data.getIndices().size();
1130 int *indices_ptr = &row_data.getIndices()[0];
1131 ierr = VecSetValues(
1133 ); CHKERRQ(ierr);
1134 } catch (const std::exception& ex) {
1135 ostringstream ss;
1136 ss << "throw in method: " << ex.what() << endl;
1137 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1138 }
1139 PetscFunctionReturn(0);
1140 }
1141 };
1142
1143 /** \brief Assemble internal force vector
1144 \ingroup gel
1145
1146 \f[
1147 (\mathbf{f}^\textrm{internal}_x)_i =
1148 \int_V
1149 \frac{\partial N_i}{\partial X_j} \sigma_{ij}
1150 \textrm{d}V
1151 \f]
1152
1153 */
1157 AssembleVector(common_data.spatialPositionName),
1158 commonData(common_data) {
1159 }
1160 PetscErrorCode doWork(
1161 int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
1162 ) {
1163 PetscFunctionBegin;
1164 try {
1165 int nb_dofs = row_data.getIndices().size();
1166 if(!nb_dofs) {
1167 PetscFunctionReturn(0);
1168 }
1169 nF.resize(nb_dofs,false);
1170 nF.clear();
1171 int nb_gauss_pts = row_data.getN().size1();
1172 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
1173 const MatrixAdaptor &diffN = row_data.getDiffN(gg,nb_dofs/3);
1174 const MatrixDouble& stress = commonData.stressTotal[gg];
1175 double val = getVolume()*getGaussPts()(3,gg);
1176 for(int dd = 0;dd<nb_dofs/3;dd++) {
1177 for(int rr = 0;rr<3;rr++) {
1178 for(int nn = 0;nn<3;nn++) {
1179 nF[3*dd+rr] += val*diffN(dd,nn)*stress(rr,nn);
1180 }
1181 }
1182 }
1183 }
1184 ierr = aSemble(row_side,row_type,row_data); CHKERRQ(ierr);
1185 } catch (const std::exception& ex) {
1186 ostringstream ss;
1187 ss << "throw in method: " << ex.what() << endl;
1188 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1189 }
1190 PetscFunctionReturn(0);
1191 }
1192 };
1193
1194 /** \brief Calculate internal forces for solvent flux
1195 \ingroup gel
1196
1197 \f[
1198 (\mathbf{f}^J_{\mu})_j =
1199 \int_V
1200 \frac{\partial N_j}{\partial X_k} J_k
1201 \textrm{d}V
1202 \f]
1203
1204 */
1208 AssembleVector(common_data.muName),
1209 commonData(common_data) {
1210 }
1211 PetscErrorCode doWork(
1212 int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
1213 ) {
1214 PetscFunctionBegin;
1215 try {
1216 int nb_dofs = row_data.getIndices().size();
1217 if(!nb_dofs) {
1218 PetscFunctionReturn(0);
1219 }
1220 nF.resize(nb_dofs,false);
1221 nF.clear();
1222 int nb_gauss_pts = row_data.getN().size1();
1223 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
1224 const MatrixAdaptor &diffN = row_data.getDiffN(gg,nb_dofs);
1225 const VectorDouble& flux = commonData.solventFlux[gg];
1226 double val = getVolume()*getGaussPts()(3,gg);
1227 for(int dd = 0;dd<nb_dofs;dd++) {
1228 for(int jj = 0;jj<3;jj++) {
1229 nF[dd] += val*diffN(dd,jj)*flux(jj);
1230 }
1231 }
1232 }
1233 ierr = aSemble(row_side,row_type,row_data); CHKERRQ(ierr);
1234 } catch (const std::exception& ex) {
1235 ostringstream ss;
1236 ss << "throw in method: " << ex.what() << endl;
1237 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1238 }
1239 PetscFunctionReturn(0);
1240 }
1241 };
1242
1243 /** \brief Calculating right hand side
1244 \ingroup gel
1245
1246 \f[
1247 (\mathbf{f}^V_\mu)_j =
1248 \int_V
1249 N_j \frac{\partial V}{\partial t}
1250 \textrm{d}V
1251 \f]
1252
1253 */
1257 AssembleVector(common_data.muName),
1258 commonData(common_data) {
1259 }
1260 PetscErrorCode doWork(
1261 int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
1262 ) {
1263 PetscFunctionBegin;
1264 try {
1265 int nb_dofs = row_data.getIndices().size();
1266 if(!nb_dofs) {
1267 PetscFunctionReturn(0);
1268 }
1269 nF.resize(nb_dofs,false);
1270 nF.clear();
1271 int nb_gauss_pts = row_data.getN().size1();
1272 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
1273 const VectorAdaptor &N = row_data.getN(gg,nb_dofs);
1274 double solvent_dot_dot = commonData.solventConcentrationDot[gg];
1275 double val = getVolume()*getGaussPts()(3,gg);
1276 nF += val*N*solvent_dot_dot;
1277 }
1278 ierr = aSemble(row_side,row_type,row_data); CHKERRQ(ierr);
1279 } catch (const std::exception& ex) {
1280 ostringstream ss;
1281 ss << "throw in method: " << ex.what() << endl;
1282 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1283 }
1284 PetscFunctionReturn(0);
1285 }
1286 };
1287
1288 /** \brief Residual strain hat
1289 \ingroup gel
1290
1291 \f[
1292 (\mathbf{r}_{\hat{\varepsilon}})_j =
1293 \int_V
1294 N_j \left(
1295 \frac{\partial \hat{\varepsilon}}{\partial t}-
1296 f(\sigma^\beta)
1297 \right)
1298 \textrm{d}V
1299 \f]
1300
1301 */
1305 AssembleVector(common_data.strainHatName),
1306 commonData(common_data) {
1307 }
1308 PetscErrorCode doWork(
1309 int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
1310 ) {
1311 PetscFunctionBegin;
1312 try {
1313 int nb_dofs = row_data.getIndices().size();
1314 if(!nb_dofs) {
1315 PetscFunctionReturn(0);
1316 }
1317 nF.resize(nb_dofs,false);
1318 nF.clear();
1319 int nb_gauss_pts = row_data.getN().size1();
1320 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
1321 const VectorAdaptor &N = row_data.getN(gg,nb_dofs/6);
1322 VectorDouble9& residual_strain_hat = commonData.residualStrainHat[gg];
1323 double val = getVolume()*getGaussPts()(3,gg);
1324 for(int dd = 0;dd<nb_dofs/6;dd++) {
1325 for(int rr = 0;rr<6;rr++) {
1326 nF[dd*6+rr] += val*N[dd]*residual_strain_hat[rr];
1327 }
1328 }
1329 }
1330 ierr = aSemble(row_side,row_type,row_data); CHKERRQ(ierr);
1331 } catch (const std::exception& ex) {
1332 ostringstream ss;
1333 ss << "throw in method: " << ex.what() << endl;
1334 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1335 }
1336 PetscFunctionReturn(0);
1337 }
1338 };
1339
1341 AssembleMatrix(string row_name,string col_name):
1342 MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
1343 row_name,col_name,UserDataOperator::OPROWCOL) {
1344 }
1345
1346 MatrixDouble K,transK;
1347 PetscErrorCode aSemble(
1348 int row_side,int col_side,
1349 EntityType row_type,EntityType col_type,
1350 DataForcesAndSourcesCore::EntData &row_data,
1351 DataForcesAndSourcesCore::EntData &col_data
1352 ) {
1353 PetscFunctionBegin;
1354 try {
1355 int nb_row = row_data.getIndices().size();
1356 int nb_col = col_data.getIndices().size();
1357 int *row_indices_ptr = &row_data.getIndices()[0];
1358 int *col_indices_ptr = &col_data.getIndices()[0];
1359 // cerr << K << endl;
1360 ierr = MatSetValues(
1361 getFEMethod()->ts_B,
1362 nb_row,row_indices_ptr,
1363 nb_col,col_indices_ptr,
1365 ); CHKERRQ(ierr);
1366 if(sYmm) {
1367 // Assemble of diagonal terms
1368 if(row_side != col_side || row_type != col_type) {
1369 transK.resize(nb_col,nb_row,false);
1370 noalias(transK) = trans(K);
1371 ierr = MatSetValues(
1372 getFEMethod()->ts_B,
1373 nb_col,col_indices_ptr,
1374 nb_row,row_indices_ptr,
1376 ); CHKERRQ(ierr);
1377 }
1378 }
1379 } catch (const std::exception& ex) {
1380 ostringstream ss;
1381 ss << "throw in method: " << ex.what() << endl;
1382 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1383 }
1384 PetscFunctionReturn(0);
1385 }
1386 };
1387
1388 /** \brief Assemble matrix \f$\mathbf{K}_{xx}\f$
1389 */
1390 struct OpLhsdxdx: public AssembleMatrix {
1392 OpLhsdxdx(CommonData &common_data):
1394 common_data.spatialPositionName,
1395 common_data.spatialPositionName
1396 ),
1397 commonData(common_data) {
1398 }
1399 MatrixDouble dStress_dx;
1400 PetscErrorCode get_dStress_dx(
1401 DataForcesAndSourcesCore::EntData &col_data,int gg
1402 ) {
1403 PetscFunctionBegin;
1404 try {
1405 int nb_col = col_data.getIndices().size();
1406 dStress_dx.resize(9,nb_col,false);
1407 dStress_dx.clear();
1408 const MatrixAdaptor diffN = col_data.getDiffN(gg,nb_col/3);
1409 MatrixDouble &jac_stress = commonData.jacStressTotal[gg];
1410 // FIXME: This can be implemented more efficiently. At this stage is
1411 // efficiency bottle neck.
1412 for(int dd = 0;dd<nb_col/3;dd++) { // DoFs in column
1413 for(int jj = 0;jj<3;jj++) { // cont. DoFs in column
1414 double a = diffN(dd,jj);
1415 for(int rr = 0;rr<3;rr++) { // Loop over dsigma_ii/dX_rr
1416 for(int ii = 0;ii<9;ii++) { // ii represents components of stress tensor
1417 dStress_dx(ii,3*dd+rr) += jac_stress(ii,3*rr+jj)*a;
1418 }
1419 }
1420 }
1421 }
1422 } catch (const std::exception& ex) {
1423 ostringstream ss;
1424 ss << "throw in method: " << ex.what() << endl;
1425 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1426 }
1427 PetscFunctionReturn(0);
1428 }
1429 PetscErrorCode doWork(
1430 int row_side,int col_side,
1431 EntityType row_type,EntityType col_type,
1432 DataForcesAndSourcesCore::EntData &row_data,
1433 DataForcesAndSourcesCore::EntData &col_data
1434 ) {
1435 PetscFunctionBegin;
1436 int nb_row = row_data.getIndices().size();
1437 int nb_col = col_data.getIndices().size();
1438 if(nb_row == 0) PetscFunctionReturn(0);
1439 if(nb_col == 0) PetscFunctionReturn(0);
1440 try {
1441 K.resize(nb_row,nb_col,false);
1442 K.clear();
1443 int nb_gauss_pts = row_data.getN().size1();
1444 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
1445 ierr = get_dStress_dx(col_data,gg); CHKERRQ(ierr);
1446 double val = getVolume()*getGaussPts()(3,gg);
1447 //cerr << dStress_dx << endl;
1448 dStress_dx *= val;
1449 const MatrixAdaptor &diffN = row_data.getDiffN(gg,nb_row/3);
1450 { //integrate element stiffness matrix
1451 for(int dd1 = 0;dd1<nb_row/3;dd1++) {
1452 for(int rr1 = 0;rr1<3;rr1++) {
1453 for(int dd2 = 0;dd2<nb_col/3;dd2++) {
1454 for(int rr2 = 0;rr2<3;rr2++) {
1455 K(3*dd1+rr1,3*dd2+rr2) += (
1456 diffN(dd1,0)*dStress_dx(3*rr1+0,3*dd2+rr2)+
1457 diffN(dd1,1)*dStress_dx(3*rr1+1,3*dd2+rr2)+
1458 diffN(dd1,2)*dStress_dx(3*rr1+2,3*dd2+rr2)
1459 );
1460 }
1461 }
1462 }
1463 }
1464 }
1465 }
1466 //cerr << "G " << getMoFEMFEPtr()->get_ref_ent() << endl << K << endl;
1467 ierr = aSemble(
1468 row_side,col_side,row_type,col_type,row_data,col_data
1469 ); CHKERRQ(ierr);
1470
1471 } catch (const std::exception& ex) {
1472 ostringstream ss;
1473 ss << "throw in method: " << ex.what() << endl;
1474 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1475 }
1476 PetscFunctionReturn(0);
1477 }
1478 };
1479
1480 /** \brief Assemble matrix \f$\mathbf{K}_{x\mu}\f$
1481 */
1484 OpLhsdxdMu(CommonData &common_data):
1486 common_data.spatialPositionName,common_data.muName
1487 ),
1488 commonData(common_data) {
1489 sYmm = false;
1490 }
1491 MatrixDouble dStress_dMu;
1492 PetscErrorCode get_dStress_dx(
1493 DataForcesAndSourcesCore::EntData &col_data,int gg
1494 ) {
1495 PetscFunctionBegin;
1496 try {
1497 int nb_col = col_data.getIndices().size();
1498 dStress_dMu.resize(9,nb_col,false);
1499 dStress_dMu.clear();
1500 const VectorAdaptor N = col_data.getN(gg);
1501 MatrixDouble &jac_stress = commonData.jacStressTotal[gg];
1502 for(int dd = 0;dd<nb_col;dd++) {
1503 double a = N[dd];
1504 for(int ii = 0;ii<9;ii++) {
1505 dStress_dMu(ii,dd) += jac_stress(ii,9+6+0)*a;
1506 }
1507 }
1508 } catch (const std::exception& ex) {
1509 ostringstream ss;
1510 ss << "throw in method: " << ex.what() << endl;
1511 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1512 }
1513 PetscFunctionReturn(0);
1514 }
1515 PetscErrorCode doWork(
1516 int row_side,int col_side,
1517 EntityType row_type,EntityType col_type,
1518 DataForcesAndSourcesCore::EntData &row_data,
1519 DataForcesAndSourcesCore::EntData &col_data
1520 ) {
1521 PetscFunctionBegin;
1522 int nb_row = row_data.getIndices().size();
1523 int nb_col = col_data.getIndices().size();
1524 if(nb_row == 0) PetscFunctionReturn(0);
1525 if(nb_col == 0) PetscFunctionReturn(0);
1526 try {
1527 K.resize(nb_row,nb_col,false);
1528 K.clear();
1529 int nb_gauss_pts = row_data.getN().size1();
1530 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
1531 ierr = get_dStress_dx(col_data,gg); CHKERRQ(ierr);
1532 double val = getVolume()*getGaussPts()(3,gg);
1533 dStress_dMu *= val;
1534 const MatrixAdaptor &diffN = row_data.getDiffN(gg,nb_row/3);
1535 { //integrate element stiffness matrix
1536 for(int dd1 = 0;dd1<nb_row/3;dd1++) {
1537 for(int rr1 = 0;rr1<3;rr1++) {
1538 for(int dd2 = 0;dd2<nb_col;dd2++) {
1539 K(3*dd1+rr1,dd2) +=
1540 diffN(dd1,0)*dStress_dMu(3*rr1+0,dd2)+
1541 diffN(dd1,1)*dStress_dMu(3*rr1+1,dd2)+
1542 diffN(dd1,2)*dStress_dMu(3*rr1+2,dd2);
1543 }
1544 }
1545 }
1546 }
1547 }
1548 ierr = aSemble(
1549 row_side,col_side,row_type,col_type,row_data,col_data
1550 ); CHKERRQ(ierr);
1551 } catch (const std::exception& ex) {
1552 ostringstream ss;
1553 ss << "throw in method: " << ex.what() << endl;
1554 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1555 }
1556 PetscFunctionReturn(0);
1557 }
1558 };
1559
1560 /** \brief Assemble matrix \f$\mathbf{K}_{x\hat{\varepsilon}}\f$
1561 */
1566 common_data.spatialPositionName,common_data.strainHatName
1567 ),
1568 commonData(common_data) {
1569 sYmm = false;
1570 }
1573 DataForcesAndSourcesCore::EntData &col_data,int gg
1574 ) {
1575 PetscFunctionBegin;
1576 try {
1577 int nb_col = col_data.getIndices().size();
1578 dStress_dStrainHat.resize(9,nb_col,false);
1579 dStress_dStrainHat.clear();
1580 const VectorAdaptor N = col_data.getN(gg);
1581 MatrixDouble &jac_stress = commonData.jacStressTotal[gg];
1582 /*cerr << N << endl;
1583 cerr << jac_stress << endl;
1584 cerr << dStress_dStrainHat << endl;
1585 cerr << nb_col << endl;*/
1586 for(int dd = 0;dd<nb_col/6;dd++) { /// DoFS in column
1587 double a = N[dd];
1588 for(int ii = 0;ii<9;ii++) { // ii for elements in stress matrix
1589 for(int rr = 0;rr<6;rr++) {
1590 dStress_dStrainHat(ii,6*dd+rr) += jac_stress(ii,9+rr)*a;
1591 }
1592 }
1593 }
1594 } catch (const std::exception& ex) {
1595 ostringstream ss;
1596 ss << "throw in method: " << ex.what() << endl;
1597 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1598 }
1599 PetscFunctionReturn(0);
1600 }
1601 PetscErrorCode doWork(
1602 int row_side,int col_side,
1603 EntityType row_type,EntityType col_type,
1604 DataForcesAndSourcesCore::EntData &row_data,
1605 DataForcesAndSourcesCore::EntData &col_data
1606 ) {
1607 PetscFunctionBegin;
1608 int nb_row = row_data.getIndices().size();
1609 int nb_col = col_data.getIndices().size();
1610 if(nb_row == 0) PetscFunctionReturn(0);
1611 if(nb_col == 0) PetscFunctionReturn(0);
1612 try {
1613 K.resize(nb_row,nb_col,false);
1614 K.clear();
1615 int nb_gauss_pts = row_data.getN().size1();
1616 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
1617 ierr = get_dStress_dStrainHat(col_data,gg); CHKERRQ(ierr);
1618 double val = getVolume()*getGaussPts()(3,gg);
1619 dStress_dStrainHat *= val;
1620 const MatrixAdaptor &diffN = row_data.getDiffN(gg,nb_row/3);
1621 { //integrate element stiffness matrix
1622 for(int dd1 = 0;dd1<nb_row/3;dd1++) {
1623 for(int rr1 = 0;rr1<3;rr1++) {
1624 for(int dd2 = 0;dd2<nb_col/6;dd2++) {
1625 for(int rr2 = 0;rr2<6;rr2++) {
1626 K(3*dd1+rr1,6*dd2+rr2) +=
1627 diffN(dd1,0)*dStress_dStrainHat(3*rr1+0,6*dd2+rr2)+
1628 diffN(dd1,1)*dStress_dStrainHat(3*rr1+1,6*dd2+rr2)+
1629 diffN(dd1,2)*dStress_dStrainHat(3*rr1+2,6*dd2+rr2);
1630 }
1631 }
1632 }
1633 }
1634 }
1635 }
1636 ierr = aSemble(
1637 row_side,col_side,row_type,col_type,row_data,col_data
1638 ); CHKERRQ(ierr);
1639 } catch (const std::exception& ex) {
1640 ostringstream ss;
1641 ss << "throw in method: " << ex.what() << endl;
1642 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1643 }
1644 PetscFunctionReturn(0);
1645 }
1646 };
1647
1648 /** \brief Assemble matrix \f$\mathbf{K}_{\hat{\varepsilon}\hat{\varepsilon}}\f$
1649 */
1654 common_data.strainHatName,common_data.strainHatName
1655 ),
1656 commonData(common_data) {
1657 }
1660 DataForcesAndSourcesCore::EntData &col_data,int gg
1661 ) {
1662 PetscFunctionBegin;
1663 try {
1664 int nb_col = col_data.getIndices().size();
1665 dStrainHat_dStrainHat.resize(6,nb_col,false);
1666 dStrainHat_dStrainHat.clear();
1667 const VectorAdaptor N = col_data.getN(gg);
1668 MatrixDouble &jac_res_strain_hat = commonData.jacStrainHat[gg];
1669 for(int dd = 0;dd<nb_col/6;dd++) { /// DoFS in column
1670 double a = N[dd];
1671 for(int ii = 0;ii<6;ii++) { // ii for elements in stress matrix
1672 for(int rr = 0;rr<6;rr++) {
1673 dStrainHat_dStrainHat(ii,6*dd+rr) += jac_res_strain_hat(ii,9+rr)*a;
1674 dStrainHat_dStrainHat(ii,6*dd+rr) += jac_res_strain_hat(ii,9+6+rr)*a*getFEMethod()->ts_a; // Time dependent element
1675 }
1676 }
1677 }
1678 } catch (const std::exception& ex) {
1679 ostringstream ss;
1680 ss << "throw in method: " << ex.what() << endl;
1681 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1682 }
1683 PetscFunctionReturn(0);
1684 }
1685 PetscErrorCode doWork(
1686 int row_side,int col_side,
1687 EntityType row_type,EntityType col_type,
1688 DataForcesAndSourcesCore::EntData &row_data,
1689 DataForcesAndSourcesCore::EntData &col_data
1690 ) {
1691 PetscFunctionBegin;
1692 int nb_row = row_data.getIndices().size();
1693 int nb_col = col_data.getIndices().size();
1694 if(nb_row == 0) PetscFunctionReturn(0);
1695 if(nb_col == 0) PetscFunctionReturn(0);
1696 try {
1697 K.resize(nb_row,nb_col,false);
1698 K.clear();
1699 int nb_gauss_pts = row_data.getN().size1();
1700 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
1701 ierr = get_dStrainHat_dStrainHat(col_data,gg); CHKERRQ(ierr);
1702 double val = getVolume()*getGaussPts()(3,gg);
1703 dStrainHat_dStrainHat *= val;
1704 const VectorAdaptor &N = row_data.getN(gg);
1705 { //integrate element stiffness matrix
1706 for(int dd1 = 0;dd1<nb_row/6;dd1++) {
1707 for(int rr1 = 0;rr1<6;rr1++) {
1708 for(int dd2 = 0;dd2<nb_col;dd2++) {
1709 K(6*dd1+rr1,dd2) += N[dd1]*dStrainHat_dStrainHat(rr1,dd2);
1710 }
1711 }
1712 }
1713 }
1714 }
1715 ierr = aSemble(
1716 row_side,col_side,row_type,col_type,row_data,col_data
1717 ); CHKERRQ(ierr);
1718 } catch (const std::exception& ex) {
1719 ostringstream ss;
1720 ss << "throw in method: " << ex.what() << endl;
1721 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1722 }
1723 PetscFunctionReturn(0);
1724 }
1725 };
1726
1727 /** \brief Assemble matrix \f$\mathbf{K}_{\hat{\varepsilon}x}\f$
1728 */
1733 common_data.strainHatName,common_data.spatialPositionName
1734 ),
1735 commonData(common_data) {
1736 sYmm = false;
1737 }
1738 MatrixDouble dStrainHat_dx;
1739 PetscErrorCode get_dStrainHat_dx(
1740 DataForcesAndSourcesCore::EntData &col_data,int gg
1741 ) {
1742 PetscFunctionBegin;
1743 try {
1744 int nb_col = col_data.getIndices().size();
1745 dStrainHat_dx.resize(6,nb_col,false);
1746 dStrainHat_dx.clear();
1747 const MatrixAdaptor diffN = col_data.getDiffN(gg,nb_col/3);
1748 MatrixDouble &jac_res_strain_hat = commonData.jacStrainHat[gg];
1749 //cerr << "a\n" << diffN << endl;
1750 //cerr << "b\n" << jac_res_strain_hat << endl;
1751 for(int dd = 0;dd<nb_col/3;dd++) { // DoFs in column
1752 for(int jj = 0;jj<3;jj++) { // cont. DoFs in column
1753 double a = diffN(dd,jj);
1754 for(int ii = 0;ii<6;ii++) { // ii for elements in strain hat
1755 for(int rr2 = 0;rr2<3;rr2++) {
1756 dStrainHat_dx(ii,3*dd+rr2) += jac_res_strain_hat(ii,3*rr2+jj)*a;
1757 }
1758 }
1759 }
1760 }
1761 } catch (const std::exception& ex) {
1762 ostringstream ss;
1763 ss << "throw in method: " << ex.what() << endl;
1764 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1765 }
1766 PetscFunctionReturn(0);
1767 }
1768 PetscErrorCode doWork(
1769 int row_side,int col_side,
1770 EntityType row_type,EntityType col_type,
1771 DataForcesAndSourcesCore::EntData &row_data,
1772 DataForcesAndSourcesCore::EntData &col_data
1773 ) {
1774 PetscFunctionBegin;
1775 int nb_row = row_data.getIndices().size();
1776 int nb_col = col_data.getIndices().size();
1777 if(nb_row == 0) PetscFunctionReturn(0);
1778 if(nb_col == 0) PetscFunctionReturn(0);
1779 try {
1780 K.resize(nb_row,nb_col,false);
1781 K.clear();
1782 int nb_gauss_pts = row_data.getN().size1();
1783 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
1784 ierr = get_dStrainHat_dx(col_data,gg); CHKERRQ(ierr);
1785 double val = getVolume()*getGaussPts()(3,gg);
1786 dStrainHat_dx *= val;
1787 //cerr << dStrainHat_dx << endl;
1788 const VectorAdaptor &N = row_data.getN(gg);
1789 { //integrate element stiffness matrix
1790 for(int dd1 = 0;dd1<nb_row/6;dd1++) {
1791 for(int rr1 = 0;rr1<6;rr1++) {
1792 for(int dd2 = 0;dd2<nb_col/3;dd2++) {
1793 for(int rr2 = 0;rr2<3;rr2++) {
1794 K(6*dd1+rr1,3*dd2+rr2) += N[dd1]*dStrainHat_dx(rr1,3*dd2+rr2);
1795 }
1796 }
1797 }
1798 }
1799 }
1800 }
1801 ierr = aSemble(
1802 row_side,col_side,row_type,col_type,row_data,col_data
1803 ); CHKERRQ(ierr);
1804 } catch (const std::exception& ex) {
1805 ostringstream ss;
1806 ss << "throw in method: " << ex.what() << endl;
1807 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1808 }
1809 PetscFunctionReturn(0);
1810 }
1811 };
1812
1813 /** \brief Assemble matrix \f$\mathbf{K}_{\mu \mu}\f$
1814 */
1819 common_data.muName,common_data.muName
1820 ),
1821 commonData(common_data) {
1822 }
1823 MatrixDouble dSolventFlux_dmu;
1824 PetscErrorCode get_dSolventFlux_dmu(
1825 DataForcesAndSourcesCore::EntData &col_data,int gg
1826 ) {
1827 PetscFunctionBegin;
1828 try {
1829 int nb_col = col_data.getIndices().size();
1830 dSolventFlux_dmu.resize(3,nb_col,false);
1831 dSolventFlux_dmu.clear();
1832 const MatrixAdaptor diffN = col_data.getDiffN(gg);
1833 MatrixDouble &jac_solvent_flux = commonData.jacSolventFlux[gg];
1834 for(int dd = 0;dd<nb_col;dd++) { // DoFs in column
1835 for(int jj = 0;jj<3;jj++) {
1836 double a = diffN(dd,jj);
1837 for(int rr2 = 0;rr2<3;rr2++) {
1838 dSolventFlux_dmu(rr2,dd) += jac_solvent_flux(rr2,jj)*a;
1839 }
1840 }
1841 }
1842 } catch (const std::exception& ex) {
1843 ostringstream ss;
1844 ss << "throw in method: " << ex.what() << endl;
1845 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1846 }
1847 PetscFunctionReturn(0);
1848 }
1849 PetscErrorCode doWork(
1850 int row_side,int col_side,
1851 EntityType row_type,EntityType col_type,
1852 DataForcesAndSourcesCore::EntData &row_data,
1853 DataForcesAndSourcesCore::EntData &col_data
1854 ) {
1855 PetscFunctionBegin;
1856 int nb_row = row_data.getIndices().size();
1857 int nb_col = col_data.getIndices().size();
1858 if(nb_row == 0) PetscFunctionReturn(0);
1859 if(nb_col == 0) PetscFunctionReturn(0);
1860 try {
1861 K.resize(nb_row,nb_col,false);
1862 K.clear();
1863 int nb_gauss_pts = row_data.getN().size1();
1864 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
1865 ierr = get_dSolventFlux_dmu(col_data,gg); CHKERRQ(ierr);
1866 double val = getVolume()*getGaussPts()(3,gg);
1867 dSolventFlux_dmu *= val;
1868 //cerr << dSolventFlux_dmu << endl;
1869 const MatrixAdaptor &diffN = row_data.getDiffN(gg);
1870 { //integrate element stiffness matrix
1871 for(int dd1 = 0;dd1<nb_row;dd1++) {
1872 for(int dd2 = 0;dd2<nb_col;dd2++) {
1873 for(int jj = 0;jj<3;jj++) {
1874 K(dd1,dd2) += diffN(dd1,jj)*dSolventFlux_dmu(jj,dd2);
1875 }
1876 }
1877 }
1878 }
1879 }
1880 ierr = aSemble(
1881 row_side,col_side,row_type,col_type,row_data,col_data
1882 ); CHKERRQ(ierr);
1883 } catch (const std::exception& ex) {
1884 ostringstream ss;
1885 ss << "throw in method: " << ex.what() << endl;
1886 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1887 }
1888 PetscFunctionReturn(0);
1889 }
1890 };
1891
1892
1893 /** \brief Assemble matrix \f$\mathbf{K}_{\mu x}\f$
1894 */
1897 OpLhsdMudx(CommonData &common_data):
1899 common_data.muName,common_data.spatialPositionName
1900 ),
1901 commonData(common_data) {
1902 sYmm = false;
1903 }
1904 VectorDouble dSolventDot_dx;
1905 PetscErrorCode get_dSolventDot_dx(
1906 DataForcesAndSourcesCore::EntData &col_data,int gg
1907 ) {
1908 PetscFunctionBegin;
1909 try {
1910 int nb_col = col_data.getIndices().size();
1911 dSolventDot_dx.resize(nb_col,false);
1912 dSolventDot_dx.clear();
1913 const MatrixAdaptor diffN = col_data.getDiffN(gg);
1914 VectorDouble &jac_solvent_rate = commonData.jacSolventDot[gg];
1915 //cerr << dSolventDot_dx << endl;
1916 //cerr << jac_solvent_rate << endl;
1917 for(int dd = 0;dd<nb_col/3;dd++) { // DoFs in column
1918 for(int jj = 0;jj<3;jj++) {
1919 double a = diffN(dd,jj);
1920 for(int rr2 = 0;rr2<3;rr2++) {
1921 dSolventDot_dx[3*dd+rr2] += jac_solvent_rate(0+3*rr2+jj)*a;
1922 dSolventDot_dx[3*dd+rr2] += jac_solvent_rate(9+3*rr2+jj)*a*getFEMethod()->ts_a;
1923 }
1924 }
1925 }
1926 //cerr << dSolventDot_dx << endl;
1927 } catch (const std::exception& ex) {
1928 ostringstream ss;
1929 ss << "throw in method: " << ex.what() << endl;
1930 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1931 }
1932 PetscFunctionReturn(0);
1933 }
1934 PetscErrorCode doWork(
1935 int row_side,int col_side,
1936 EntityType row_type,EntityType col_type,
1937 DataForcesAndSourcesCore::EntData &row_data,
1938 DataForcesAndSourcesCore::EntData &col_data
1939 ) {
1940 PetscFunctionBegin;
1941 int nb_row = row_data.getIndices().size();
1942 int nb_col = col_data.getIndices().size();
1943 if(nb_row == 0) PetscFunctionReturn(0);
1944 if(nb_col == 0) PetscFunctionReturn(0);
1945 try {
1946 K.resize(nb_row,nb_col,false);
1947 K.clear();
1948 int nb_gauss_pts = row_data.getN().size1();
1949 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
1950 ierr = get_dSolventDot_dx(col_data,gg); CHKERRQ(ierr);
1951 double val = getVolume()*getGaussPts()(3,gg);
1952 dSolventDot_dx *= val;
1953 //cerr << dSolventDot_dx << endl;
1954 const VectorAdaptor &N = row_data.getN(gg);
1955 { //integrate element stiffness matrix
1956 for(int dd1 = 0;dd1<nb_row;dd1++) {
1957 for(int dd2 = 0;dd2<nb_col;dd2++) {
1958 K(dd1,dd2) += N[dd1]*dSolventDot_dx[dd2];
1959 }
1960 }
1961 }
1962 }
1963 ierr = aSemble(
1964 row_side,col_side,row_type,col_type,row_data,col_data
1965 ); CHKERRQ(ierr);
1966 } catch (const std::exception& ex) {
1967 ostringstream ss;
1968 ss << "throw in method: " << ex.what() << endl;
1969 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1970 }
1971 PetscFunctionReturn(0);
1972 }
1973 };
1974
1975 /**
1976 * \brief Used to post proc stresses, energy, source term.
1977 *
1978 * calculation of principal stresses. for e.g.
1979 *
1980 */
1982
1984 moab::Interface &postProcMesh;
1985 std::vector<EntityHandle> &mapGaussPts;
1986
1988 const string field_name,
1989 moab::Interface &post_proc_mesh,
1990 std::vector<EntityHandle> &map_gauss_pts,
1991 CommonData &common_data
1992 ):
1993 MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(field_name,UserDataOperator::OPROW),
1994 commonData(common_data),
1995 postProcMesh(post_proc_mesh),
1996 mapGaussPts(map_gauss_pts) {
1997 // Set that this operator is executed only for vertices on element
1998 doVertices = true;
1999 doEdges = false;
2001 doTris = false;
2002 doTets = false;
2003 doPrisms = false;
2004 }
2005
2006 PetscErrorCode doWork(
2007 int side,EntityType type,DataForcesAndSourcesCore::EntData &data
2008 ) {
2009 //
2010
2011 PetscFunctionBegin;
2012
2013 if(type!=MBVERTEX) PetscFunctionReturn(9);
2014 double def_VAL[9];
2015 bzero(def_VAL,9*sizeof(double));
2016 Tag th_stress_total;
2017 rval = postProcMesh.tag_get_handle(
2018 "STRESS_TOTAL",9,MB_TYPE_DOUBLE,th_stress_total,MB_TAG_CREAT|MB_TAG_SPARSE,def_VAL
2019 ); CHKERRQ_MOAB(rval);
2020 Tag th_flux;
2021 rval = postProcMesh.tag_get_handle(
2023 ); CHKERRQ_MOAB(rval);
2024
2025 const int nb_gauss_pts = mapGaussPts.size();
2026 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
2027 rval = postProcMesh.tag_set_data(th_stress_total,&mapGaussPts[gg],1,&(commonData.stressTotal[gg](0,0))); CHKERRQ_MOAB(rval);
2028 rval = postProcMesh.tag_set_data(th_flux,&mapGaussPts[gg],1,&(commonData.solventFlux[gg][0])); CHKERRQ_MOAB(rval);
2029 }
2030
2031 PetscFunctionReturn(0);
2032
2033 }
2034
2035 };
2036
2037
2039
2041 string pRoblem;
2042 string fE;
2045 vector<int> tAgs;
2046
2048
2049 bool iNit;
2050 int pRT;
2051
2053 MoFEM::Interface &m_field,
2054 string problem,
2055 string fe,
2056 CommonData &common_data,
2058 vector<int> tags
2059 ):
2060 FEMethod(),
2061 mField(m_field),
2062 pRoblem(problem),
2063 fE(fe),
2064 commonData(common_data),
2065 cE(ce),
2066 tAgs(tags),
2067 postProc(m_field),
2068 iNit(false) {
2069
2070 PetscBool flg = PETSC_TRUE;
2071 ierr = PetscOptionsGetInt(
2072 PETSC_NULL,PETSC_NULL,"-my_output_prt",&pRT,&flg
2073 ); CHKERRABORT(PETSC_COMM_WORLD,ierr);
2074 if(flg!=PETSC_TRUE) {
2075 pRT = 1;
2076 }
2077 }
2078
2079 PetscErrorCode preProcess() {
2080 PetscFunctionBegin;
2081 PetscFunctionReturn(0);
2082 }
2083
2084 PetscErrorCode operator()() {
2085 PetscFunctionBegin;
2086 PetscFunctionReturn(0);
2087 }
2088
2089 PetscErrorCode postProcess() {
2090 PetscFunctionBegin;
2091
2092 if(!iNit) {
2093
2099 postProc.getOpPtrVector().push_back
2100 (new Gel::OpGetDataAtGaussPts("SPATIAL_POSITION",commonData,false,true));
2101 postProc.getOpPtrVector().push_back
2102 (new Gel::OpGetDataAtGaussPts("SPATIAL_POSITION_DOT",commonData,false,true));
2103 postProc.getOpPtrVector().push_back
2105 postProc.getOpPtrVector().push_back
2106 (new Gel::OpGetDataAtGaussPts("HAT_EPS",commonData,true,false,MBTET));
2107 postProc.getOpPtrVector().push_back
2108 (new Gel::OpGetDataAtGaussPts("HAT_EPS_DOT",commonData,true,false,MBTET));
2109 postProc.getOpPtrVector().push_back
2110 (new Gel::OpJacobian("SPATIAL_POSITION",tAgs,cE,commonData,true,false));
2111 postProc.getOpPtrVector().push_back
2113
2114
2115 iNit = true;
2116 }
2117 int step;
2118 ierr = TSGetTimeStepNumber(ts,&step); CHKERRQ(ierr);
2119
2120 if((step)%pRT==0) {
2122 ostringstream sss;
2123 sss << "out_" << step << ".h5m";
2124 rval = postProc.postProcMesh.write_file(sss.str().c_str(),"MOAB","PARALLEL=WRITE_PART"); CHKERRQ_MOAB(rval);
2125 }
2126 PetscFunctionReturn(0);
2127 }
2128
2129 };
2130
2131};
2132
2133}
2134
2135#endif //__GEL_HPP__
2136
2137/***************************************************************************//**
2138 * \defgroup gel Gel model
2139 * \ingroup user_modules
2140 ******************************************************************************/
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
static PetscErrorCode ierr
@ COL
Definition: definitions.h:123
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:454
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
constexpr int order
@ F
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
constexpr auto field_name
const int N
Definition: speed_test.cpp:3
AssembleMatrix(string row_name, string col_name)
Definition: Gels.hpp:1341
PetscErrorCode aSemble(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: Gels.hpp:1347
AssembleVector(string field_name)
Definition: Gels.hpp:1119
PetscErrorCode aSemble(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: Gels.hpp:1124
Gel material parameters.
Definition: Gels.hpp:88
double vIscosity
Viscosity of Gel.
Definition: Gels.hpp:98
double gAlpha
Shear modulus spring alpha.
Definition: Gels.hpp:92
double gBeta
Shear modulus spring beta.
Definition: Gels.hpp:94
double vBeta
Poisson ratio spring beta.
Definition: Gels.hpp:93
double gBetaHat
Shear modulus dashpot.
Definition: Gels.hpp:96
double pErmeability
Permeability of Gel.
Definition: Gels.hpp:97
double oMega
Volume per solvent molecule.
Definition: Gels.hpp:99
double vBetaHat
Poisson ratio dashpot.
Definition: Gels.hpp:95
double mU0
Definition: Gels.hpp:100
Range tEts
Test in the block.
Definition: Gels.hpp:102
double vAlpha
Poisson ratio spring alpha.
Definition: Gels.hpp:91
Common data for gel model.
Definition: Gels.hpp:345
vector< VectorDouble > jacSolventDot
Definition: Gels.hpp:365
vector< MatrixDouble > jacSolventFlux
Definition: Gels.hpp:364
vector< MatrixDouble3by3 > stressTotal
Definition: Gels.hpp:357
map< int, int > nbActiveResults
Definition: Gels.hpp:369
vector< VectorDouble3 > solventFlux
Definition: Gels.hpp:358
map< string, vector< VectorDouble > > dataAtGaussPts
Definition: Gels.hpp:354
vector< double * > jacRowPtr
Definition: Gels.hpp:362
map< int, int > nbActiveVariables
Definition: Gels.hpp:369
map< string, vector< MatrixDouble > > gradAtGaussPts
Definition: Gels.hpp:355
vector< MatrixDouble > jacStressTotal
Definition: Gels.hpp:363
vector< MatrixDouble > jacStrainHat
Definition: Gels.hpp:366
vector< double > solventConcentrationDot
Definition: Gels.hpp:359
vector< VectorDouble9 > residualStrainHat
Definition: Gels.hpp:360
Constitutive model functions.
Definition: Gels.hpp:114
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > stressBetaHat
Stress as result of volume change due to solvent concentration.
Definition: Gels.hpp:142
ublas::vector< TYPE, ublas::bounded_array< TYPE, 3 > > gradientMu
Definition: Gels.hpp:131
virtual PetscErrorCode calculateResidualStrainHat()
Definition: Gels.hpp:301
virtual PetscErrorCode calculateStressBeta()
Calculate stress in spring beta.
Definition: Gels.hpp:226
ConstitutiveEquation(map< int, BlockMaterialData > &data)
Definition: Gels.hpp:119
virtual PetscErrorCode calculateSolventFlux()
Calculate flux.
Definition: Gels.hpp:319
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > residualStrainHat
Residual for calculation epsilon hat.
Definition: Gels.hpp:154
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > FDot
Definition: Gels.hpp:127
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > strainTotal
Total strain applied at integration point.
Definition: Gels.hpp:137
virtual PetscErrorCode calculateStrainHatFlux()
Calculate rate of strain hat.
Definition: Gels.hpp:255
virtual PetscErrorCode calculateSolventConcentrationDot()
Calculate solvent concentration rate.
Definition: Gels.hpp:332
virtual PetscErrorCode calculateStrainTotal()
Calculate total strain.
Definition: Gels.hpp:166
ublas::vector< TYPE, ublas::bounded_array< TYPE, 9 > > solventFlux
Solvent flux.
Definition: Gels.hpp:152
TYPE solventConcentrationDot
Volume rate change.
Definition: Gels.hpp:153
virtual PetscErrorCode calculateCauchyDefromationTensor()
Definition: Gels.hpp:156
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > stressTotal
Total stress.
Definition: Gels.hpp:151
map< int, BlockMaterialData > & dAta
Definition: Gels.hpp:117
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > stressAlpha
Stress generated by spring alpha.
Definition: Gels.hpp:139
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > strainHatDot
Internal variable, strain in dashpot beta.
Definition: Gels.hpp:129
virtual PetscErrorCode calculateStressBetaHat()
Calculate stress due to concentration of solvent molecules.
Definition: Gels.hpp:280
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > strainHatFlux
Rate of dashpot (beta) strain.
Definition: Gels.hpp:141
virtual PetscErrorCode calculateStressTotal()
Definition: Gels.hpp:292
virtual PetscErrorCode calculateStressAlpha()
Calculate stress in spring alpha.
Definition: Gels.hpp:200
TYPE mU
Solvent concentration.
Definition: Gels.hpp:130
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > strainHat
Internal variable, strain in dashpot beta.
Definition: Gels.hpp:128
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > F
Definition: Gels.hpp:126
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > gradientU
Definition: Gels.hpp:136
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > C
Cauchy deformation.
Definition: Gels.hpp:135
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > stressBeta
Stress generated by spring beta.
Definition: Gels.hpp:140
virtual PetscErrorCode calculateTraceStrainTotalDot()
Definition: Gels.hpp:181
definition of volume element
Definition: Gels.hpp:379
PetscErrorCode preProcess()
function is run at the beginning of loop
Definition: Gels.hpp:394
CommonData & commonData
Definition: Gels.hpp:381
int getRule(int order)
Definition: Gels.hpp:390
Takes into account HO geometry.
Definition: Gels.hpp:382
GelFE(MoFEM::Interface &m_field, CommonData &common_data)
Definition: Gels.hpp:384
PetscErrorCode postProcess()
function is run at the end of loop
Definition: Gels.hpp:416
Definition: Gels.hpp:2038
PetscErrorCode postProcess()
function is run at the end of loop
Definition: Gels.hpp:2089
PetscErrorCode operator()()
function is run for every finite element
Definition: Gels.hpp:2084
vector< int > tAgs
Definition: Gels.hpp:2045
string fE
Definition: Gels.hpp:2042
bool iNit
Definition: Gels.hpp:2049
boost::shared_ptr< Gel::ConstitutiveEquation< adouble > > & cE
Definition: Gels.hpp:2044
PostProcVolumeOnRefinedMesh postProc
Definition: Gels.hpp:2047
MoFEM::Interface & mField
Definition: Gels.hpp:2040
MonitorPostProc(MoFEM::Interface &m_field, string problem, string fe, CommonData &common_data, boost::shared_ptr< Gel::ConstitutiveEquation< adouble > > &ce, vector< int > tags)
Definition: Gels.hpp:2052
int pRT
Definition: Gels.hpp:2050
PetscErrorCode preProcess()
function is run at the beginning of loop
Definition: Gels.hpp:2079
string pRoblem
Definition: Gels.hpp:2041
CommonData & commonData
Definition: Gels.hpp:2043
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator field value.
Definition: Gels.hpp:470
OpGetDataAtGaussPts(const string field_name, CommonData &common_data, bool calc_val, bool calc_grad, EntityType zero_at_type=MBVERTEX)
Definition: Gels.hpp:451
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: Gels.hpp:1076
PetscErrorCode calculateAtIntPtsSolventFlux()
Definition: Gels.hpp:908
PetscErrorCode calculateAtIntPtsStressTotal()
Definition: Gels.hpp:833
PetscErrorCode recordStressTotal()
Definition: Gels.hpp:577
VectorDouble activeVariables
Definition: Gels.hpp:575
boost::shared_ptr< Gel::ConstitutiveEquation< adouble > > cE
Definition: Gels.hpp:544
PetscErrorCode calculateJacobian(TagEvaluate te)
Definition: Gels.hpp:808
vector< int > tagS
Definition: Gels.hpp:543
PetscErrorCode calculateAtIntPtsSolventDot()
Definition: Gels.hpp:958
PetscErrorCode calculateAtIntPtrsResidualStrainHat()
Definition: Gels.hpp:1006
CommonData & commonData
Definition: Gels.hpp:545
PetscErrorCode recordSolventConcentrationDot()
Definition: Gels.hpp:677
PetscErrorCode recordSolventFlux()
Definition: Gels.hpp:641
map< int, int > & nbActiveResults
Definition: Gels.hpp:551
PetscErrorCode recordResidualStrainHat()
Definition: Gels.hpp:719
map< int, int > & nbActiveVariables
Definition: Gels.hpp:550
OpJacobian(const string field_name, vector< int > tags, boost::shared_ptr< Gel::ConstitutiveEquation< adouble > > ce, CommonData &common_data, bool calculate_residual, bool calculate_jacobian)
Definition: Gels.hpp:553
PetscErrorCode calculateFunction(TagEvaluate te, double *ptr)
Definition: Gels.hpp:786
Assemble matrix .
Definition: Gels.hpp:1815
MatrixDouble dSolventFlux_dmu
Definition: Gels.hpp:1823
OpLhsdMudMu(CommonData &common_data)
Definition: Gels.hpp:1817
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: Gels.hpp:1849
PetscErrorCode get_dSolventFlux_dmu(DataForcesAndSourcesCore::EntData &col_data, int gg)
Definition: Gels.hpp:1824
CommonData & commonData
Definition: Gels.hpp:1816
Assemble matrix .
Definition: Gels.hpp:1895
VectorDouble dSolventDot_dx
Definition: Gels.hpp:1904
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: Gels.hpp:1934
CommonData & commonData
Definition: Gels.hpp:1896
OpLhsdMudx(CommonData &common_data)
Definition: Gels.hpp:1897
PetscErrorCode get_dSolventDot_dx(DataForcesAndSourcesCore::EntData &col_data, int gg)
Definition: Gels.hpp:1905
OpLhsdStrainHatdStrainHat(CommonData &common_data)
Definition: Gels.hpp:1652
PetscErrorCode get_dStrainHat_dStrainHat(DataForcesAndSourcesCore::EntData &col_data, int gg)
Definition: Gels.hpp:1659
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: Gels.hpp:1685
OpLhsdStrainHatdx(CommonData &common_data)
Definition: Gels.hpp:1731
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: Gels.hpp:1768
PetscErrorCode get_dStrainHat_dx(DataForcesAndSourcesCore::EntData &col_data, int gg)
Definition: Gels.hpp:1739
Assemble matrix .
Definition: Gels.hpp:1482
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: Gels.hpp:1515
OpLhsdxdMu(CommonData &common_data)
Definition: Gels.hpp:1484
MatrixDouble dStress_dMu
Definition: Gels.hpp:1491
CommonData & commonData
Definition: Gels.hpp:1483
PetscErrorCode get_dStress_dx(DataForcesAndSourcesCore::EntData &col_data, int gg)
Definition: Gels.hpp:1492
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: Gels.hpp:1601
OpLhsdxdStrainHat(CommonData &common_data)
Definition: Gels.hpp:1564
PetscErrorCode get_dStress_dStrainHat(DataForcesAndSourcesCore::EntData &col_data, int gg)
Definition: Gels.hpp:1572
Assemble matrix .
Definition: Gels.hpp:1390
OpLhsdxdx(CommonData &common_data)
Definition: Gels.hpp:1392
CommonData & commonData
Definition: Gels.hpp:1391
MatrixDouble dStress_dx
Definition: Gels.hpp:1399
PetscErrorCode get_dStress_dx(DataForcesAndSourcesCore::EntData &col_data, int gg)
Definition: Gels.hpp:1400
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: Gels.hpp:1429
Used to post proc stresses, energy, source term.
Definition: Gels.hpp:1981
std::vector< EntityHandle > & mapGaussPts
Definition: Gels.hpp:1985
CommonData & commonData
Definition: Gels.hpp:1983
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Gels.hpp:2006
moab::Interface & postProcMesh
Definition: Gels.hpp:1984
OpPostProcGel(const string field_name, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, CommonData &common_data)
Definition: Gels.hpp:1987
Calculating right hand side.
Definition: Gels.hpp:1254
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: Gels.hpp:1260
OpRhsSolventConcetrationDot(CommonData &common_data)
Definition: Gels.hpp:1256
Calculate internal forces for solvent flux.
Definition: Gels.hpp:1205
OpRhsSolventFlux(CommonData &common_data)
Definition: Gels.hpp:1207
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: Gels.hpp:1211
Residual strain hat.
Definition: Gels.hpp:1302
OpRhsStrainHat(CommonData &common_data)
Definition: Gels.hpp:1304
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: Gels.hpp:1308
Assemble internal force vector.
Definition: Gels.hpp:1154
OpRhsStressTotal(CommonData &common_data)
Definition: Gels.hpp:1156
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: Gels.hpp:1160
Implementation of Gel constitutive model.
Definition: Gels.hpp:74
GelFE feLhs
Definition: Gels.hpp:436
Gel(MoFEM::Interface &m_field)
Definition: Gels.hpp:438
GelFE feRhs
Definition: Gels.hpp:436
map< int, BlockMaterialData > blockMaterialData
Definition: Gels.hpp:105
MoFEM::Interface & mFiled
Definition: Gels.hpp:83
@ RESIDUALSTRAINHAT
Definition: Gels.hpp:80
CommonData commonData
Definition: Gels.hpp:376
boost::shared_ptr< Gel::ConstitutiveEquation< adouble > > constitutiveEquationPtr
Definition: Gels.hpp:340
const Problem * problemPtr
raw pointer to problem
bool & doTris
\deprectaed
bool & doPrisms
\deprectaed
bool sYmm
If true assume that matrix is symmetric structure.
bool & doVertices
\deprectaed If false skip vertices
bool & doTets
\deprectaed
bool & doEdges
\deprectaed If false skip edges
\deprectaed
Deprecated interface functions.
structure for User Loop Methods on finite elements
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
@ OPCOL
operator doWork function is executed on FE columns
@ 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
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
virtual MoFEMErrorCode postProcess()
function is run at the end of loop
TS ts
time solver
Vec & ts_F
residual vector
TSContext ts_ctx
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
Vec & ts_u_t
time derivative of state vector
Mat & ts_B
Preconditioner for ts_A.
VolumeElementForcesAndSourcesCore * getVolumeFE() const
return pointer to Generic Volume Finite Element object
VolumeElementForcesAndSourcesCore(Interface &m_field, const EntityType type=MBTET)
std::vector< EntityHandle > mapGaussPts
moab::Interface & postProcMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Post processing.