NeoPZ
swelling.cpp
Go to the documentation of this file.
1 
6 #include "swelling.h"
7 #include "pzelmat.h"
8 #include "pzbndcond.h"
9 #include "pzmatrix.h"
10 #include "pzfmatrix.h"
11 #include "pzerror.h"
12 #include "checkconv.h"
13 #include <math.h>
14 
15 using namespace std;
16 
17 #include <cmath>
18 
19 #ifdef _AUTODIFF
20 REAL TPZSwelling::gFaraday = 96.4853;
21 REAL TPZSwelling::gVPlus = 2.3;
22 REAL TPZSwelling::gVMinus = 15.17;
23 REAL TPZSwelling::gRGas = 8.3145;
24 REAL TPZSwelling::gTemp = 293.;
25 REAL TPZSwelling::gMuRef[3] = {0.,0.,0.};
26 #else
27 STATE TPZSwelling::gFaraday = 96.4853;
28 STATE TPZSwelling::gVPlus = 2.3;
29 STATE TPZSwelling::gVMinus = 15.17;
30 STATE TPZSwelling::gRGas = 8.3145;
31 STATE TPZSwelling::gTemp = 293.;
32 STATE TPZSwelling::gMuRef[3] = {0.,0.,0.};
33 #endif
34 
35 TPZFMatrix<STATE> TPZSwelling::gState;
36 TPZFMatrix<REAL> TPZSwelling::gphi(1,1,1.);
37 TPZFMatrix<REAL> TPZSwelling::gdphi(3,1,0.34);
38 
39 #ifdef _AUTODIFF
40 
41 void ToMatrix(TPZVec<FADREAL> &vec, TPZFMatrix<STATE> &ek);
42 
43 #endif
44 
45 TPZSwelling::TPZSwelling(int matindex, STATE lambda, STATE shear, STATE alfa, STATE M, STATE Gamma, STATE Kperm, STATE DPlus, STATE DMinus,
46  STATE rHinder, STATE Cfc, STATE Nf0, STATE NPlus0, STATE NMinus0) :
47 TPZRegisterClassId(&TPZSwelling::ClassId),
48 TPZMaterial(matindex) {
49  fComputationMode = 0;
50  fLambda = lambda;
51  fShear = shear;
52  fAlfa = alfa;
53  fM = M;
54  fGamma = Gamma;
55  fKperm.Redim(3,3);
56  fDPlus = DPlus;
57  fDMinus = DMinus;
58  frHinder = rHinder;
59  fInitDeform = 1.;
60  fCfc = Cfc;
61  fNf0 = Nf0;
62  fNPlus0 = NPlus0;
63  fNMinus0 = NMinus0;
64  int ic;
65  for(ic=0; ic<3; ic++) {
66  fKperm(ic,ic) = Kperm;
67  }
68  fTheta = 1.;
69  fDelt = 0.5;
70 }
71 
72 TPZSwelling::~TPZSwelling() {
73 }
74 
75 
76 void TPZSwelling::Print(std::ostream &out) {
77  TPZMaterial::Print(out);
78  out << "name of material : " << Name() << "\n";
79  out << "properties : \n";
80  out << "Computation mode : " << fComputationMode << endl;
81  out << "Compression modulus " << fLambda << endl;
82  out << "Shear modulus " << fShear << endl;
83  out << "Biot coupling coeficient " << fAlfa << endl;
84  out << "Storage modulus " << fM << endl;
85  out << "Osmotic coeficient " << fGamma << endl;
86  out << "Hydraulic permeability " << fKperm << endl;
87  out << "Diffusion coeficient for cations " << fDPlus << endl;
88  out << "Diffusion coeficient for anions " << fDMinus << endl;
89  out << "Hindrance factor " << frHinder << endl;
90  out << "Initial deformation " << fInitDeform << endl;
91  out << "Fixed charge density " << fCfc << endl;
92  out << "Initial fluid volume fraction " << fNf0 << endl;
93  out << "Initial cation volume fraction " << fNPlus0 << endl;
94  out << "Initial anion volume fraction " << fNMinus0 << endl;
95  out << "Weight factor for time integration of ionic conservation law " << fTheta << endl;
96  out << "Timestep " << fDelt << endl;
97  out << "Faraday constant " << gFaraday << endl;
98  out << "Molar volume cation " <<gVPlus << endl;
99  out << "Molar volume anions " << gVMinus << endl;
100  out << "Gas constant " << gRGas << endl;
101  out << "Absolute temperature " << gTemp << endl;
102  int ieq;
103  for(ieq=0; ieq<3; ieq++) {
104  out << "Reference chemical potentials [" << ieq << "] " << gMuRef[ieq] << endl;
105  }
106 }
107 
108 
109 
111 int TPZSwelling::VariableIndex(const std::string &name){
112  return TPZMaterial::VariableIndex(name);
113 }
114 
115 int TPZSwelling::NSolutionVariables(int var){
116 
117  return TPZMaterial::NSolutionVariables(var);
118 }
119 
120 void TPZSwelling::Solution(TPZVec<STATE> &Sol,TPZFMatrix<STATE> &DSol,TPZFMatrix<REAL> &axes,int var,TPZVec<STATE> &Solout){
121 
122  TPZMaterial::Solution(Sol,DSol,axes,var,Solout);
123 }
124 
125 
126 
127 
128 void TPZSwelling::ContributeBC(TPZMaterialData &data,
129  REAL weight,
130  TPZFMatrix<STATE> &ek,
131  TPZFMatrix<STATE> &ef,
132  TPZBndCond &bc) {
133 
134  TPZFMatrix<REAL> &phi = data.phi;
135  int numbersol = data.sol.size();
136  if (numbersol != 1) {
137  DebugStop();
138  }
139 
140  TPZVec<STATE> &sol=data.sol[0];
141 
142  if(bc.Material() != this){
143  PZError << "TPZMatHyperElastic.ContributeBC : this material don't exists \n";
144  }
145 
146  if(bc.Type() < 0 && bc.Type() > 2){
147  PZError << "ContributeBC.aplybc, unknown boundary condition type : "<<bc.Type() << endl;
148  }
149 
150  int ndof = NStateVariables();
151  int nnod = ek.Rows()/ndof;
152  int r = ndof;
153 
154  int idf,jdf,in,jn;
155  switch(bc.Type()){
156  case 0:
157  for(in=0 ; in<nnod ; ++in){
158  for(idf = 0;idf<r;idf++) {
159  (ef)(in*r+idf,0) += gBigNumber*phi(in,0)*(bc.Val2()(idf,0)-sol[idf])*weight;
160  }
161  for(jn=0 ; jn<nnod ; ++jn) {
162  for(idf = 0;idf<r;idf++) {
163  ek(in*r+idf,jn*r+idf) += gBigNumber*phi(in,0)*phi(jn,0)*weight;
164  }
165  }
166  }
167  break;
168 
169  case 1:
170  for(in=0 ; in<nnod ; ++in){
171  for(idf = 0;idf<r;idf++) {
172  //(ef)(in*r+idf,0) += weight*phi(in,0)*(bc.Val2()(idf,0)-sol[idf]);
173  (ef)(in*r+idf,0) += weight*phi(in,0)*(bc.Val2()(idf,0));
174  }
175  }
176  break;
177 
178  case 2:
179  for(in=0 ; in<nnod ; ++in){
180  for(idf = 0;idf<r;idf++) {
181  for (jdf=0; jdf<r; jdf++){
182  (ef)(in*r+idf,0) += phi(in,0)*bc.Val1()(idf,jdf)*(bc.Val2()(jdf,0)-sol[jdf])*weight;
183  }
184  for(jn=0 ; jn<nnod ; ++jn) {
185  for(idf = 0;idf<r;idf++) {
186  for(jdf = 0;jdf<r;jdf++) {
187  ek(in*r+idf,jn*r+jdf) += bc.Val1()(idf,jdf)*phi(in,0)*phi(jn,0)*weight;
188  }
189  }
190  }
191  }
192 
193  }
194 
195  }//fim switch
196 }
197 
198 #ifdef _AUTODIFF
199 
200 /* The function below makes the correspondence between the dsol vector
201  and a matrix ordered F operator
202  */
203 inline int ith(const int i, const int j)
204 {
205  return i*3+j;
206 }
207 
208 
209 void TPZSwelling::ContributeElastEnergy(TPZVec<FADFADREAL> &dsol,
210  FADFADREAL &U, REAL weight)
211 {
212  FADFADREAL J, TrC; // J = det(F); TrC = Trace(C)
213  int numeq = dsol[0].size();
214  FADREAL defaultFAD(numeq, REAL(0.), REAL(0.));
215  FADFADREAL defaultFADFAD(numeq, defaultFAD, defaultFAD);
216  TPZManVector<FADFADREAL,9> GradMap(9);
217  int iel;
218  for(iel=0; iel<9; iel++) {
219  GradMap[iel] = dsol[iel]*FADREAL(fInitDeform);
220  }
221 
222  for(iel=0; iel<3; iel++) GradMap[ith(iel,iel)].val().val() += fInitDeform;
223 
224  TrC = GradMap[0]*GradMap[0]+GradMap[1]*GradMap[1]+GradMap[2]*GradMap[2]+
225  GradMap[3]*GradMap[3]+GradMap[4]*GradMap[4]+GradMap[5]*GradMap[5]+
226  GradMap[6]*GradMap[6]+GradMap[7]*GradMap[7]+GradMap[8]*GradMap[8];
227 
228  DebugStop();
229  /*
230  J = GradMap[ith(0,0)] * GradMap[ith(1,1)] * GradMap[ith(2,2)] +
231  GradMap[ith(0,1)] * GradMap[ith(1,2)] * GradMap[ith(2,0)] +
232  GradMap[ith(0,2)] * GradMap[ith(1,0)] * GradMap[ith(2,1)] -
233  GradMap[ith(0,2)] * GradMap[ith(1,1)] * GradMap[ith(2,0)] -
234  GradMap[ith(0,1)] * GradMap[ith(1,0)] * GradMap[ith(2,2)] -
235  GradMap[ith(0,0)] * GradMap[ith(1,2)] * GradMap[ith(2,1)]; // J = det(F)
236 
237  // expression of the elastic energy
238  U += (J*J - FADREAL(1.)) * FADREAL(weight*fLambda/4.) -
239  log( J ) * FADREAL(weight*(fLambda/2.+fShear)) +
240  (TrC - FADREAL(3.)) * FADREAL(weight*fShear/2.);
241  */
242 }
243 
247 void TPZSwelling::ContributeResidual(TPZVec<REAL> & x,
248  TPZVec<FADREAL> & sol,
249  TPZVec<FADREAL> &dsol,
250  TPZFMatrix<REAL> &phi,
251  TPZFMatrix<REAL> &dphi,
252  TPZVec<FADREAL> &RES,
253  REAL weight){
254  const int nstate = 8;
255 
256  if(fComputationMode == 0) {
257  ContributePrevResidual(x,sol,dsol,phi,dphi,RES,weight);
258  return;
259  }
260  int numeq = sol[0].size();
261  FADREAL defaultFAD(numeq, REAL(0.), REAL(0.));
262  FADFADREAL defaultFADFAD(numeq, defaultFAD, defaultFAD);
263  FADFADREAL deform(defaultFADFAD);
264  TPZManVector<FADFADREAL> dsolFADFAD(9,defaultFADFAD);
265 
266  // extend the gradient of the solution to a variable which contains second derivatives
267  int ieq,der;
268  for(der=0; der<9; der++) {
269  dsolFADFAD[der].val().val() = dsol[der].val();
270  for(ieq=0; ieq<numeq; ieq++) {
271  dsolFADFAD[der].val().fastAccessDx(ieq) = dsol[der].fastAccessDx(ieq);
272  dsolFADFAD[der].fastAccessDx(ieq).val() = dsol[der].fastAccessDx(ieq);
273  }
274  }
275  // The first derivative of the deform variable contains the residual vector
276  // The second derivative contains the tangent matrix
277  ContributeElastEnergy(dsolFADFAD, deform, weight);
278 
279  int jeq;
280  for(ieq=0; ieq<numeq; ieq++) {
281  RES[ieq].val() += deform.fastAccessDx(ieq).val();
282  for(jeq=0; jeq<numeq; jeq++) {
283  RES[ieq].fastAccessDx(jeq) += deform.fastAccessDx(ieq).fastAccessDx(jeq);
284  }
285  }
286 
287  //return;
288  // cout << "dsol " << dsol;
289 
290  // compute the gradient of the map induced by the displacement of the element
291  // include the derivative of the map with respect to the solution
292  FADREAL GradMap[3][3];
293  GradMap[0][0] = dsol[0]+REAL(1.);
294  GradMap[0][1] = dsol[1];
295  GradMap[0][2] = dsol[2];
296  GradMap[1][0] = dsol[3];
297  GradMap[1][1] = dsol[4]+REAL(1.);
298  GradMap[1][2] = dsol[5];
299  GradMap[2][0] = dsol[6];
300  GradMap[2][1] = dsol[7];
301  GradMap[2][2] = dsol[8]+REAL(1.);
302 
303  DebugStop();
304  /*
305  // compute the determinant of the map
306  // this computation will carry the derivative w.r.t the solution
307  FADREAL J = GradMap[0][0] * GradMap[1][1] * GradMap[2][2] +
308  GradMap[0][1] * GradMap[1][2] * GradMap[2][0] +
309  GradMap[0][2] * GradMap[1][0] * GradMap[2][1] -
310  GradMap[0][2] * GradMap[1][1] * GradMap[2][0] -
311  GradMap[0][1] * GradMap[1][0] * GradMap[2][2] -
312  GradMap[0][0] * GradMap[1][2] * GradMap[2][1]; // J = det(F)
313 
314  // compute the inverse of the map, including its derivatives
315  FADREAL GradMapInv[3][3];
316  for(ieq=0; ieq<3; ieq++) {
317  int ieqp = (ieq+1)%3;
318  int ieqpp = (ieq+2)%3;
319  for(jeq=0; jeq<3; jeq++) {
320  int jeqp = (jeq+1)%3;
321  int jeqpp = (jeq+2)%3;
322  GradMapInv[ieq][jeq] = (GradMap[ieqp][jeqp]*GradMap[ieqpp][jeqpp]-GradMap[ieqpp][jeqp]*GradMap[ieqp][jeqpp])/J;
323  // cout << ieq << ' ' << jeq << endl << ieqp << ' ' << jeqp << '+' << ieqpp << ' ' << jeqpp << endl
324  // << ieqpp << ' ' << jeqp << '-' << ieqp << ' ' << jeqpp << endl << GradMapInv[ieq][jeq];
325  }
326  }
327 
328  // compute the Lagrangian volume fractions in function of mu, pres and ksi
329  TPZVec<FADREAL> N(3);
330  ComputeN(sol,N);
331 
332  int ishape,nshape = phi.Rows();
333  for(ishape=0; ishape<nshape; ishape++) {
334  // compute the gradient of the weighting function in the deformed configuration
335  FADREAL GradPhi[3];
336  for(ieq=0; ieq<3; ieq++) {
337  GradPhi[ieq] = GradMapInv[0][ieq]*dphi(0,ishape)+GradMapInv[1][ieq]*dphi(1,ishape)+GradMapInv[2][ieq]*dphi(2,ishape);
338  // Add the contribution of the pressure to the linear momentum equation
339  RES[ishape*nstate+ieq] += sol[3]*(GradPhi[ieq]*weight)*J;
340  }
341  // Add the contribution of the rate of change of the volume to the mass balance of the mixture equation
342  RES[ishape*nstate+3]+= J*(phi(ishape,0)*weight);
343  // Add the contribution of the Lagrange volume fractions to the mass balance of the mixture equation
344  RES[ishape*nstate+3] -= (N[0]+N[1]+N[2])*(phi(ishape,0)*weight);
345  // Add the contribution to electro neutrality equation
346  RES[ishape*nstate+7] -= ((N[1]/(FADREAL)gVPlus)-(N[2]/(FADREAL)gVMinus))*(phi(ishape,0)*gFaraday*weight);
347  for(ieq=0; ieq<3; ieq++) {
348  // Add the contribution to the mass balance of the constituents
349  RES[ishape*nstate+4+ieq] += (N[ieq]*phi(ishape,0)*weight);
350  STATE KGradPhi = 0.;
351  for(jeq=0; jeq<3; jeq++) {
352  KGradPhi += dphi(jeq,ishape)*fKperm(ieq,jeq)*fTheta*fDelt*weight;
353  }
354  // Add the contribution of the difusion of the constituents
355  RES[ishape*nstate+4+ieq] += (FADREAL)KGradPhi*dsol[ieq+3*(4+ieq)];
356  }
357  }
358  */
359 }
360 
361 void TPZSwelling::ContributePrevResidual(TPZVec<REAL> & x,
362  TPZVec<FADREAL> & sol,
363  TPZVec<FADREAL> &dsol,
364  TPZFMatrix<REAL> &phi,
365  TPZFMatrix<REAL> &dphi,
366  TPZVec<FADREAL> &RES,
367  REAL weight){
368  const int nstate = 8;
369 
370  if(fComputationMode != 0) {
371  cout << "TPZSwelling::ContributePrevResidual should not be called\n";
372  return;
373  }
374  int numeq = sol[0].size();
375  FADREAL defaultFAD(numeq, REAL(0.), REAL(0.));
376  FADFADREAL defaultFADFAD(numeq, defaultFAD, defaultFAD);
377  FADFADREAL deform(defaultFADFAD);
378  TPZManVector<FADFADREAL> dsolFADFAD(9,defaultFADFAD);
379 
380  // extend the gradient of the solution to a variable which contains second derivatives
381  int ieq,der;
382  for(der=0; der<9; der++) {
383  dsolFADFAD[der].val().val() = dsol[der].val();
384  for(ieq=0; ieq<numeq; ieq++) {
385  dsolFADFAD[der].val().fastAccessDx(ieq) = dsol[der].fastAccessDx(ieq);
386  dsolFADFAD[der].fastAccessDx(ieq).val() = dsol[der].fastAccessDx(ieq);
387  }
388  }
389  // The first derivative of the deform variable contains the residual vector
390  // The second derivative contains the tangent matrix
391 
392  int jeq;
393 
394  // compute the gradient of the map induced by the displacement of the element
395  // include the derivative of the map with respect to the solution
396  REAL GradMap[3][3];
397  GradMap[0][0] = dsol[0].val()+1.;
398  GradMap[0][1] = dsol[1].val();
399  GradMap[0][2] = dsol[2].val();
400  GradMap[1][0] = dsol[3].val();
401  GradMap[1][1] = dsol[4].val()+1.;
402  GradMap[1][2] = dsol[5].val();
403  GradMap[2][0] = dsol[6].val();
404  GradMap[2][1] = dsol[7].val();
405  GradMap[2][2] = dsol[8].val()+1.;
406 
407  DebugStop();
408  /*
409  // compute the determinant of the map
410  // this computation will carry the derivative w.r.t the solution
411  REAL J = GradMap[0][0] * GradMap[1][1] * GradMap[2][2] +
412  GradMap[0][1] * GradMap[1][2] * GradMap[2][0] +
413  GradMap[0][2] * GradMap[1][0] * GradMap[2][1] -
414  GradMap[0][2] * GradMap[1][1] * GradMap[2][0] -
415  GradMap[0][1] * GradMap[1][0] * GradMap[2][2] -
416  GradMap[0][0] * GradMap[1][2] * GradMap[2][1]; // J = det(F)
417 
418  // compute the Lagrangian volume fractions in function of mu, pres and ksi
419  TPZVec<REAL> N(3);
420  ComputeN(sol,N);
421 
422  int ishape,nshape = phi.Rows();
423  for(ishape=0; ishape<nshape; ishape++) {
424  // Add the contribution of the rate of change of the volume to the mass balance of the mixture equation
425  RES[ishape*nstate+3] -= J*(phi(ishape,0)*weight);
426  // Add the contribution of the Lagrange volume fractions to the mass balance of the mixture equation
427  RES[ishape*nstate+3] += (N[0]+N[1]+N[2])*(phi(ishape,0)*weight);
428  // Add the contribution to electro neutrality equation
429  RES[ishape*nstate+7] += (N[1]/gVPlus-N[2]/gVMinus)*(phi(ishape,0)*gFaraday*weight);
430  for(ieq=0; ieq<3; ieq++) {
431  // Add the contribution to the mass balance of the constituents
432  RES[ishape*nstate+4+ieq] += (N[ieq]*phi(ishape,0)*weight);
433  REAL KGradPhi = 0.;
434  for(jeq=0; jeq<3; jeq++) {
435  KGradPhi += dphi(jeq,ishape)*fKperm(ieq,jeq)*(1.-fTheta)*fDelt*weight;
436  }
437  // Add the contribution of the difusion of the constituents
438  RES[ishape*nstate+4+ieq] += KGradPhi*dsol[ieq+3*(4+ieq)];
439  }
440  }
441  */
442 }
443 
444 void TPZSwelling::ComputeW(FADFADREAL &W, TPZVec<STATE> &N) {
445  DebugStop();
446  /*
447  FADREAL defaultFAD(3,REAL(0.),REAL(0.));
448  FADFADREAL defaultFADFAD(3,defaultFAD,defaultFAD);
449  W = defaultFADFAD;
450  FADFADREAL NFAD[3] = {defaultFADFAD,defaultFADFAD,defaultFADFAD};
451 
452  // transform the N (volume fractions) in second order derivatives
453  int ieq;
454  for(ieq=0; ieq<3; ieq++) {
455  NFAD[ieq].val().val() = N[ieq];
456  NFAD[ieq].val().fastAccessDx(ieq) = 1.;
457  NFAD[ieq].fastAccessDx(ieq).val() = 1.;
458  }
459  W += FADREAL(gMuRef[0])*NFAD[0]+FADREAL(gMuRef[1])*NFAD[1]+
460  FADREAL(gMuRef[2])*NFAD[2]-
461  FADREAL(gRGas*gTemp*fGamma)*(NFAD[1]*FADREAL(1./gVPlus)+NFAD[2]*FADREAL(1./gVMinus))*log(NFAD[0])+
462  FADREAL(gRGas*gTemp/gVPlus)*NFAD[1]*(log(NFAD[1]/FADREAL(gVPlus))-FADREAL(1.))+
463  FADREAL(gRGas*gTemp/gVMinus)*NFAD[2]*(log(NFAD[2]/FADREAL(gVMinus))-FADREAL(1.));
464  */
465 }
466 
467 void TPZSwelling::ComputeN(TPZVec<STATE> &mu, STATE ksi, STATE pressure, TPZVec<STATE> &N) {
468 
469  TPZFMatrix<STATE> res,tangent;
470  ExactSolution(mu,ksi,pressure,N);
471  NResidual(mu,ksi,pressure,N,res,tangent);
472  // compute the residual of the current N values
473  REAL resnorm = Norm(res);
474  int ieq;
475  while(resnorm > 1.e-6) {
476  // Compute the correction as suggested by the Newton method
477  // This only works if the residual is suficiently small
478  tangent.SolveDirect(res,ELU);
479  for(ieq=0; ieq<3; ieq++) {
480  N[ieq] -= res(ieq,0);
481  }
482  NResidual(mu,ksi,pressure,N,res,tangent);
483  resnorm = Norm(res);
484  }
485 }
486 
487 void TPZSwelling::NResidual(TPZVec<STATE> &mu, STATE ksi, STATE pressure, TPZVec<STATE> &N, TPZFMatrix<STATE> &res, TPZFMatrix<STATE> &tangent) {
488 
489  FADREAL defaultFAD(3,REAL(0.),STATE(0.));
490  FADFADREAL defaultFADFAD(3,defaultFAD,defaultFAD),W;
491  W = defaultFADFAD;
492 
493  // Compute the mixing energy as a second order derivative variable
494  // the gradient contributes to the residual
495  // the second derivatives contribute to the tangent matrix
496  ComputeW(W,N);
497 
498  res.Redim(3,1);
499  tangent.Redim(3,3);
500  STATE z[3] = {0.,1.,-1.};
501  STATE v[3] = {1.,gVPlus,gVMinus};
502  int ieq,jeq;
503  for(ieq=0; ieq<3; ieq++) {
504  res(ieq,0) = -mu[ieq]+W.val().fastAccessDx(ieq)+gFaraday*z[ieq]*ksi/v[ieq]+pressure;
505  for(jeq=0; jeq<3; jeq++) {
506  tangent(ieq,jeq) = W.fastAccessDx(ieq).fastAccessDx(jeq);
507  }
508  }
509 }
510 
511 void TPZSwelling::NResidual(TPZVec<FADREAL> &sol, TPZVec<FADREAL> &N) {
512 
513  // compute the values of N numerically and its derivatives by applying a single newton iteration using FAD variables
514  TPZManVector<STATE,3> muloc(3), Nloc(3);
515  // create the residual vector and tangent matrix for solving for N
516  TPZVec<FADREAL> res(3);
517  TPZVec<TPZVec<FADREAL> > tangent(3,res);
518 
519  FADFADREAL W;
520 
521  int ieq;
522  REAL z[3] = {0.,1.,-1.};
523  REAL v[3] = {1.,gVPlus,gVMinus};
524  FADREAL &ksi = sol[7];
525  FADREAL &pressure = sol[3];
526  for(ieq=0; ieq<3; ieq++) {
527  muloc[ieq]= sol[ieq+4].val();
528  Nloc[ieq] = N[ieq].val();
529  }
530  // compute the value of N by any method (without carrying the derivatives)
531  ExactSolution(muloc,ksi.val(),pressure.val(),Nloc);
532  ComputeW(W,Nloc);
533 
534  // initialize the N variable with the correct solution but with zero derivatives
535  for(ieq=0; ieq<3; ieq++) {
536  N[ieq] = FADREAL(sol[0].size(),((REAL)Nloc[ieq]),0.);
537  }
538  // perform a single Newton iteration
539  // build the tangent matrix and residual (including the derivatives)
540  int jeq;
541  for(ieq=0; ieq<3; ieq++) {
542  res[ieq] = -sol[ieq+4]+W.val().fastAccessDx(ieq)+gFaraday*z[ieq]*ksi/v[ieq]+pressure;
543  for(jeq=0; jeq<3; jeq++) {
544  tangent[ieq][jeq] = FADREAL(sol[0].size(),W.fastAccessDx(ieq).fastAccessDx(jeq),0);
545  }
546  }
547  // update by a single Newton iteration
548  Solve(tangent,res);
549  for(ieq=0; ieq<3; ieq++) {
550  N[ieq] -= res[ieq];
551  }
552 }
553 
554 void TPZSwelling::Solve(TPZVec<TPZVec<FADREAL> > &tangent, TPZVec<FADREAL> &res) {
555  // implements a simple LU decomposition using FADREAL variables
556  int neq = res.NElements();
557  int k,i,j;
558  for(k=0; k<neq; k++) {
559  FADREAL pivot = tangent[k][k];
560  for ( i = k+1; i < neq; i++ ) {
561  tangent[i][k] /= pivot;
562  for (j = k+1; j < neq; j++ ) tangent[i][j] -= tangent[i][k]*tangent[k][j];
563  res[i] -= tangent[i][k]*res[k];
564  }
565  }
566  for ( i = neq-1; i >= 0; i-- ) {
567  for (j = i+1; j < neq ; j++ ) {
568  res[i] -= tangent[i][j]*res[j];
569  }
570  res[i] /= tangent[i][i];
571  }
572 }
573 
574 int TPZSwelling::main() {
575 
576  // program created for testing purposes
577  int matindex = 1;
578  STATE lambda = 1;
579  STATE shear = 0.5;
580  STATE alfa = 1.;
581  STATE M = 0.2;
582  STATE Gamma = 0.9;
583  STATE Kperm = 1.e-4;
584  STATE DPlus = 5.e-4;
585  STATE DMinus = 5.e-4;
586  STATE rHinder = 0.4;
587  STATE Cfc = -2.e-4;
588  STATE Nf0 = 0.8;
589  // initiele externe zoutconcentratie
590  STATE C0 = 0.15e-3;
591  STATE cplus0 = (-Cfc+sqrt(Cfc*Cfc+4.*C0*C0))/2.;
592  STATE cminus0 = (Cfc+sqrt(Cfc*Cfc+4.*C0*C0))/2.;
593  STATE NPlus0 = gVPlus*cplus0*Nf0;
594  STATE NMinus0 = gVMinus*cminus0*Nf0;
595 
596  TPZSwelling test(matindex, lambda, shear, alfa, M, Gamma, Kperm, DPlus, DMinus,
597  rHinder, Cfc, Nf0, NPlus0, NMinus0);
598 
599  test.Print(cout);
600 
601  // initialize a set of variables as they would be setup by the element
602  int ieq,d;
603  TPZFMatrix<REAL> phi(1,1),dphi(3,1);
604  phi(0,0) = 1.;
605  for(ieq=0; ieq<3; ieq++) {
606  dphi(ieq,0) = 0.34;
607  }
608  FADREAL defaultFAD(8,REAL(0.),REAL(0.));
609  TPZVec<FADREAL> sol(8,defaultFAD),dsol(24,defaultFAD);
610  TPZVec<REAL> x(3,1.);
611  REAL weight = 0.33;
612  TPZVec<FADREAL> RES(8,defaultFAD);
613  REAL values[8] = {1.,1.,1.,0.001,-0.7243,-9258.,-1401.,-15};
614 
615  TPZVec<STATE> mu(3),N(3,0.);
616  for(ieq=0; ieq<3; ieq++) mu[ieq] = values[ieq+4];
617  STATE ksi,J = 1.0062,pres;
618  test.ComputeInitialGuess(mu,J,pres,ksi,N);
619  values[3] = pres;
620  values[7] = ksi;
621  cout << "Based on J " << J << " and mu \n" << mu << "\nThe initial guesses pres " << pres << " ksi " << ksi << "\n N\n" << N << endl;
622 
623 
624  for(ieq=0; ieq<8; ieq++) {
625  sol[ieq].val() = values[ieq];
626  sol[ieq].fastAccessDx(ieq) = phi(0);
627  for(d=0; d<3; d++) {
628  dsol[ieq*3 + d].val() = dphi(d,0)*values[ieq];
629  dsol[ieq*3 + d].fastAccessDx(ieq) = dphi(d,0);
630  }
631  }
632 
633  // procedure to check whether the stiffness matrix is tangent to the residual vector
634 
635  STATE rangeval[11] = {0.1,0.1,0.1,0.0001,0.01,1.,1.,0.1,0.,0.,0.};
636  TPZFMatrix<STATE> state(11,1),range(11,1,0.0);
637  for(ieq=0; ieq<8; ieq++) {
638  state(ieq,0) = values[ieq];
639  range(ieq,0) = rangeval[ieq];
640  }
641  for(ieq=8; ieq<11; ieq++) {
642  state(ieq,0) = N[ieq-8];
643  range(ieq,0) = rangeval[ieq];
644  }
645  TPZVec<REAL> coefs(1,1.);
646  CheckConvergence(test,state,range,coefs);
647 
648  // sample computation of a contribution to a stiffness matrix
649  // this procedure was built to check whether the stiffness matrix is symetric
650 
651  test.ContributeResidual(x,sol,dsol,phi,dphi,RES,weight);
652  TPZFMatrix<STATE> ek;
653  ToMatrix(RES,ek);
654  ek.Print("stiffness matrix");
655 
656  cout << RES;
657  return 0;
658 
659 }
660 
661 void ToMatrix(TPZVec<FADREAL> &vec, TPZFMatrix<STATE> &ek) {
662  int nel = vec.NElements();
663  ek.Redim(nel,nel);
664  int ieq,jeq;
665  for(ieq=0; ieq<nel; ieq++) {
666  for(jeq=0; jeq<nel; jeq++) {
667  ek(ieq,jeq) = vec[ieq].fastAccessDx(jeq);
668  }
669  }
670 }
671 
672 /* Methods to implement automatic conformity checks */
673 int TPZSwelling::NumCases() {
674  return 1;
675 }
676 
677 void TPZSwelling::LoadState(TPZFMatrix<STATE> &state) {
678  if(state.Rows() != 11) {
679  cout << "TPZSwelling LoadState wrong number of variables expect bad things\n";
680  }
681  gState = state;
682 }
683 void TPZSwelling::ComputeTangent(TPZFMatrix<STATE> &tangent,TPZVec<REAL> &coefs, int cases) {
684  tangent.Redim(11,11);
685  int ieq,jeq,d;
686  FADREAL defaultFAD(8,REAL(0.),REAL(0.));
687  TPZVec<FADREAL> sol(8,defaultFAD),dsol(24,defaultFAD);
688  TPZVec<REAL> x(3,1.);
689  REAL weight = 0.33;
690  TPZVec<FADREAL> RES(8,defaultFAD);
691  for(ieq=0; ieq<8; ieq++) {
692  sol[ieq].val() = gphi(0,0)*gState(ieq,0);
693  sol[ieq].fastAccessDx(ieq) = gphi(0);
694  for(d=0; d<3; d++) {
695  dsol[ieq*3 + d].val() = gdphi(d,0)*gState(ieq,0);
696  dsol[ieq*3 + d].fastAccessDx(ieq) = gdphi(d,0);
697  }
698  }
699  ContributeResidual(x,sol,dsol,gphi,gdphi,RES,weight);
700  for(ieq=0; ieq<8; ieq++) {
701  for(jeq=0; jeq<8; jeq++) {
702  tangent(ieq,jeq) = RES[ieq].fastAccessDx(jeq);
703  }
704  }
705  TPZVec<STATE> N(3,0.);
706  for(ieq=0; ieq<3; ieq++) N[ieq] = gState(8+ieq,0);
707  FADREAL defwFAD(3,REAL(0.),REAL(0.));
708  FADFADREAL defwFADFAD(3,defwFAD,defwFAD);
709  FADFADREAL W(defwFADFAD);
710  ComputeW(W,N);
711  for(ieq=0; ieq<3; ieq++) {
712  for(jeq=0; jeq<3; jeq++) {
713  tangent(8+ieq,8+jeq) = W.fastAccessDx(ieq).fastAccessDx(jeq);
714  }
715  }
716 }
717 void TPZSwelling::Residual(TPZFMatrix<STATE> &res, int cases) {
718  res.Redim(11,1);
719  int ieq,d;
720  FADREAL defaultFAD(8,REAL(0.),REAL(0.));
721  TPZVec<FADREAL> sol(8,defaultFAD),dsol(24,defaultFAD);
722  TPZVec<REAL> x(3,1.);
723  REAL weight = 0.33;
724  TPZVec<FADREAL> RES(8,defaultFAD);
725  for(ieq=0; ieq<8; ieq++) {
726  sol[ieq].val() = gphi(0,0)*gState(ieq,0);
727  sol[ieq].fastAccessDx(ieq) = gphi(0);
728  for(d=0; d<3; d++) {
729  dsol[ieq*3 + d].val() = gdphi(d,0)*gState(ieq,0);
730  dsol[ieq*3 + d].fastAccessDx(ieq) = gdphi(d,0);
731  }
732  }
733  ContributeResidual(x,sol,dsol,gphi,gdphi,RES,weight);
734  for(ieq=0; ieq<8; ieq++) {
735  res(ieq,0) = RES[ieq].val();
736  }
737  TPZVec<STATE> N(3,0.);
738  for(ieq=0; ieq<3; ieq++) N[ieq] = gState(8+ieq,0);
739  FADREAL defwFAD(3,REAL(0.),REAL(0.));
740  FADFADREAL defwFADFAD(3,defwFAD,defwFAD);
741  FADFADREAL W(defwFADFAD);
742  ComputeW(W,N);
743  for(ieq=0; ieq<3; ieq++) {
744  res(8+ieq,0) = W.fastAccessDx(ieq).val();
745  }
746 }
747 
748 #endif
749 
750 void TPZSwelling::ComputeInitialGuess(TPZVec<STATE> &mu, STATE J, STATE &pres, STATE &ksi, TPZVec<STATE> &N) {
751 
752  pres = 0.;
753  STATE expcontr = exp((gVPlus*(mu[1]-gMuRef[1]-pres)+gVMinus*(mu[2]-gMuRef[2]-pres))/(gRGas*gTemp));
754  STATE expcontr2 = expcontr*STATE(4.);
755  STATE sqrtval = sqrt(fCfc*fCfc+expcontr2);
756  STATE cplus = (-fCfc + sqrtval)/2.;
757  STATE cmin = (fCfc+sqrtval)/2.;
758 
759  int iter;
760  for(iter=0; iter<5; iter++) {
761  pres = (mu[0]-gMuRef[0])+fGamma*gRGas*gTemp*(cplus+cmin);
762  cout.precision(12);
763  cout << "pres " << pres << " cplus " << cplus << " cmin " << cmin << " expcontr " << expcontr << endl;
764  STATE expcontr = exp((gVPlus*(mu[1]-gMuRef[1]-pres)+gVMinus*(mu[2]-gMuRef[2]-pres))/(gRGas*gTemp));
765  STATE expcontr2 = expcontr*STATE(4.);
766  sqrtval = sqrt(fCfc*fCfc+expcontr2);
767  cplus = (-fCfc + sqrtval)/2.;
768  cmin = (fCfc+sqrtval)/2.;
769  }
770  cout.precision(12);
771  cout << "cplus " << cplus << " cmin " << cmin << endl;
772  ksi= (gRGas*gTemp*log(cmin/cplus)+gVPlus*(mu[1]-gMuRef[1]-pres)-gVMinus*(mu[2]-gMuRef[2]-pres))/(2.*gFaraday);
773  N[0] = J-1+fNf0;
774  N[1] = N[0]*gVPlus*cplus;
775  N[2] = N[0]*gVMinus*cmin;
776 
777 }
778 
779 void TPZSwelling::ExactSolution(TPZVec<STATE> &mu, STATE ksi, STATE pres, TPZVec<STATE> &N) {
780  STATE c1,c2;
781  c1 = (mu[1]-gFaraday*ksi/gVPlus - pres)*gVPlus/(gRGas*gTemp);
782  c2 = (mu[2]+gFaraday*ksi/gVMinus - pres)*gVMinus/(gRGas*gTemp);
783 
784  STATE fac = exp(c1)+exp(c2);
785  STATE fac2 = (-mu[0]+pres)/(gRGas*gTemp*fGamma);
786  N[0] = pow(fac2/fac,1./(fGamma-1));
787 
788  N[1] = pow(N[0],fGamma)*exp(c1)*gVPlus;
789  N[2] = pow(N[0],fGamma)*exp(c2)*gVMinus;
790 
791 }
792 
793 #ifdef _AUTODIFF
794 
795 void TPZSwelling::ComputeN(TPZVec<FADREAL> &sol, TPZVec<FADREAL> &N) {
796 
797  // computes the analytic solution of N carrying the derivatives
798  FADREAL &pres = sol[3];
799  FADREAL &mu0 = sol[4];
800  FADREAL &mu1 = sol[5];
801  FADREAL &mu2 = sol[6];
802  FADREAL &ksi = sol[7];
803 
804  FADREAL expc1,expc2;
805  expc1 = exp((mu1-gFaraday*ksi/gVPlus - pres)*gVPlus/(gRGas*gTemp));
806 
807  expc2 = exp((mu2+gFaraday*ksi/gVMinus - pres)*gVMinus/(gRGas*gTemp));
808 
809  FADREAL fac2 = (-mu0+pres)/(gRGas*gTemp*fGamma);
810  N[0] = pow(fac2/(expc1+expc2),1./(fGamma-1));
811  FADREAL N0Gamma = pow(N[0].val(),(REAL)fGamma);
812 
813  N[1] = N0Gamma*expc1*gVPlus;
814  N[2] = N0Gamma*expc2*gVMinus;
815 }
816 
817 void TPZSwelling::ComputeN(TPZVec<FADREAL> &sol, TPZVec<REAL> &N) {
818 
819  // computes the analytic solution of N carrying the derivatives
820 
821  REAL &pres = sol[3].val();
822  REAL &mu0 = sol[4].val();
823  REAL &mu1 = sol[5].val();
824  REAL &mu2 = sol[6].val();
825  REAL &ksi = sol[7].val();
826 
827  REAL expc1,expc2;
828  expc1 = exp((mu1-gFaraday*ksi/gVPlus - pres)*gVPlus/(gRGas*gTemp));
829 
830  expc2 = exp((mu2+gFaraday*ksi/gVMinus - pres)*gVMinus/(gRGas*gTemp));
831 
832  REAL fac2 = (-mu0+pres)/(gRGas*gTemp*fGamma);
833  N[0] = pow(fac2/(expc1+expc2),((REAL)(1./(fGamma-1))));
834  REAL N0Gamma = pow(N[0],((REAL)fGamma));
835 
836  N[1] = N0Gamma*expc1*gVPlus;
837  N[2] = N0Gamma*expc2*gVMinus;
838 }
839 
840 #endif
841 
842 int TPZSwelling::ClassId() const{
843  return Hash("TPZSwelling") ^ TPZMaterial::ClassId() << 1;
844 }
TPZSwelling::fDPlus
STATE fDPlus
Diffusion coeficient for cations [mm^2/s].
Definition: swelling.h:42
TPZSwelling::~TPZSwelling
virtual ~TPZSwelling()
Definition: swelling.cpp:72
TPZMaterial::Solution
virtual void Solution(TPZMaterialData &data, int var, TPZVec< STATE > &Solout)
Returns the solution associated with the var index based on the finite element approximation.
Definition: TPZMaterial.cpp:200
TPZSwelling::fKperm
TPZFMatrix< STATE > fKperm
Hydraulic permeability .
Definition: swelling.h:30
CheckConvergence
void CheckConvergence(TConv &obj, TPZFMatrix< STATE > &state, TPZFMatrix< STATE > &range, TPZVec< REAL > &coefs)
Implements a general procedure to check whether the class TConv implements a consistente tangent matr...
Definition: checkconv.h:23
TPZSwelling::fNMinus0
STATE fNMinus0
Initial anion volume fraction (no dimension)
Definition: swelling.h:56
TPZSwelling::Solution
virtual void Solution(TPZVec< STATE > &Sol, TPZFMatrix< STATE > &DSol, TPZFMatrix< REAL > &axes, int var, TPZVec< STATE > &Solout) override
Computes a post-processed solution variable corresponding to the variable index.
Definition: swelling.cpp:120
TPZManVector
Implements a vector class which allows to use external storage provided by the user. Utility.
Definition: pzquad.h:16
TPZSwelling::fM
STATE fM
Storage modulus [N/mm^2].
Definition: swelling.h:40
TPZSwelling::fDMinus
STATE fDMinus
Diffusion coeficient for anions [mm^2/s].
Definition: swelling.h:44
bc
clarg::argBool bc("-bc", "binary checkpoints", false)
TPZSwelling::NumCases
int NumCases()
Methods needed to perform convergence checks.
TPZSwelling::Print
virtual void Print(std::ostream &out) override
Prints out the data associated with the material.
Definition: swelling.cpp:76
TPZBndCond::Type
int Type()
Definition: pzbndcond.h:249
TPZSwelling::fShear
STATE fShear
Shear modulus [N/mm^2].
Definition: swelling.h:36
test
Definition: test.py:1
TPZMaterial::VariableIndex
virtual int VariableIndex(const std::string &name)
Returns the variable index associated with the name.
Definition: TPZMaterial.cpp:141
pzerror.h
Defines PZError.
TPZSwelling
The TPZSwelling class implements a numerical model of swelling material coupling flow through porous ...
Definition: swelling.h:21
TPZSwelling::fDelt
STATE fDelt
Timestep [s].
Definition: swelling.h:61
TPZSwelling::NStateVariables
virtual int NStateVariables() const override
Number of state variables, in this case: 3 displacements, 1 pressure, 3 eletrochemical potencials...
Definition: swelling.h:127
TPZVec
This class implements a simple vector storage scheme for a templated class T. Utility.
Definition: pzgeopoint.h:19
val
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
TPZSwelling::gdphi
static TPZFMatrix< REAL > gdphi
Definition: swelling.h:281
TPZBndCond::Val2
TPZFMatrix< STATE > & Val2(int loadcase=0)
Definition: pzbndcond.h:255
TPZMaterialData::phi
TPZFNMatrix< 220, REAL > phi
vector of shapefunctions (format is dependent on the value of shapetype)
Definition: pzmaterialdata.h:55
TPZSwelling::ExactSolution
void ExactSolution(TPZVec< STATE > &mu, STATE ksi, STATE pres, TPZVec< STATE > &N)
Definition: swelling.cpp:779
TPZSwelling::ComputeTangent
void ComputeTangent(TPZFMatrix< STATE > &tangent, TPZVec< REAL > &coefs, int cases)
Computes the tangent matrix for a given loadcase.
TPZSwelling::fAlfa
STATE fAlfa
Biot coupling coeficient (no dimension)
Definition: swelling.h:38
TPZMaterial::Print
virtual void Print(std::ostream &out=std::cout)
Prints out the data associated with the material.
Definition: TPZMaterial.cpp:103
TPZMaterial
This abstract class defines the behaviour which each derived class needs to implement.
Definition: TPZMaterial.h:39
pzelmat.h
Contains declaration of TPZElementMatrix struct which associates an element matrix with the coeficien...
TPZSwelling::gTemp
static STATE gTemp
Absolute temperature [K].
Definition: swelling.h:89
pzbndcond.h
Contains the TPZBndCond class which implements a boundary condition for TPZMaterial objects...
TPZSwelling::gRGas
static STATE gRGas
gas constant [Nmm/(mmol K)]
Definition: swelling.h:87
TPZSwelling::fCfc
STATE fCfc
Fixed charge density [mmoleq/mm^3].
Definition: swelling.h:50
TPZVec::size
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
TPZSwelling::fInitDeform
STATE fInitDeform
Initial deformation (assuming everything occurs isotropically, a constant is suficient (no dimension)...
Definition: swelling.h:48
TPZSwelling::fNf0
STATE fNf0
Initial fluid volume fraction (do dimension)
Definition: swelling.h:52
Norm
TVar Norm(const TPZFMatrix< TVar > &A)
Returns the norm of the matrix A.
pzfmatrix.h
Contains TPZMatrixclass which implements full matrix (using column major representation).
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
TPZBndCond
This class defines the boundary condition for TPZMaterial objects.
Definition: pzbndcond.h:29
TPZSwelling::Name
virtual std::string Name() override
Returns the name of the material.
Definition: swelling.h:131
ELU
Definition: pzmatrix.h:52
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
TPZSwelling::fTheta
STATE fTheta
Timestepping parameter theta (no dimension)
Definition: swelling.h:59
sqrt
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ sqrt
Definition: tfadfunc.h:120
TPZSwelling::gFaraday
static STATE gFaraday
Faraday constant [C/mmol].
Definition: swelling.h:81
TPZSwelling::ComputeN
void ComputeN(TPZVec< STATE > &N, TPZVec< STATE > &mu, STATE ksi)
Computes the value of the N coeficients in function of ksi and mus, iterative method, inverting the Hessian of W.
test.res
string res
Definition: test.py:151
TPZSwelling::fComputationMode
int fComputationMode
Computation mode: 0->residual wrt to time 1->residual and tangent wrt to time .
Definition: swelling.h:28
TPZSwelling::gphi
static TPZFMatrix< REAL > gphi
Variables which holds the state variables used in the check convergence procedure.
Definition: swelling.h:281
TPZBndCond::Val1
TPZFMatrix< STATE > & Val1()
Definition: pzbndcond.h:253
TPZSwelling::fGamma
STATE fGamma
Osmotic coeficient (no dimension)
Definition: swelling.h:34
TPZFMatrix::Redim
int Redim(const int64_t newRows, const int64_t newCols) override
Redimension a matrix and ZERO your elements.
Definition: pzfmatrix.h:616
TPZSwelling::TPZSwelling
TPZSwelling(int matindex, STATE lambda, STATE shear, STATE alfa, STATE M, STATE Gamma, STATE Kperm, STATE DPlus, STATE DMinus, STATE rHinder, STATE Cfc, STATE Nf0, STATE NPlus0, STATE NMinus0)
Constructor of the class, where the user needs to specify the most important parameters.
Definition: swelling.cpp:45
TPZMaterial::gBigNumber
static REAL gBigNumber
Big number to penalization method, used for Dirichlet conditions.
Definition: TPZMaterial.h:83
TPZSwelling::LoadState
void LoadState(TPZFMatrix< STATE > &state)
Loads the state within the current object, to be used when computing the tangent matrix.
Hash
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
pzmatrix.h
Contains TPZMatrix<TVar>class, root matrix class.
TPZMaterial::ClassId
int ClassId() const override
Unique identifier for serialization purposes.
Definition: TPZMaterial.cpp:382
TPZSwelling::fNPlus0
STATE fNPlus0
Initial cation volume fraction (no dimension)
Definition: swelling.h:54
TPZSwelling::ComputeInitialGuess
void ComputeInitialGuess(TPZVec< STATE > &mu, STATE J, STATE &pres, STATE &ksi, TPZVec< STATE > &N)
Computes the aproximate values of the pressure, ksi and N based on mu and J by direct inversion of th...
Definition: swelling.cpp:750
pow
TPZFlopCounter pow(const TPZFlopCounter &orig, const TPZFlopCounter &xp)
Returns the power and increments the counter of the power.
Definition: pzreal.h:487
idf
static int idf[6]
Definition: TPZCompMeshTools.cpp:135
log
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ log
Definition: tfadfunc.h:130
TPZMaterial::NSolutionVariables
virtual int NSolutionVariables(int var)
Returns the number of variables associated with the variable indexed by var.
Definition: TPZMaterial.cpp:176
TPZSwelling::VariableIndex
virtual int VariableIndex(const std::string &name) override
Definition: swelling.cpp:111
TPZSwelling::gMuRef
static STATE gMuRef[3]
Reference chemical potentials (order f,plus,minus) [mV].
Definition: swelling.h:91
TPZSwelling::gVPlus
static STATE gVPlus
Molar volume cation [mm^3/mmol].
Definition: swelling.h:83
TPZSwelling::ClassId
virtual int ClassId() const override
Define the class id associated with the class.
Definition: swelling.cpp:842
TPZSwelling::frHinder
STATE frHinder
Hindrance factor (no dimension)
Definition: swelling.h:46
exp
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ exp
Definition: tfadfunc.h:125
TPZMatrix::Print
virtual void Print(std::ostream &out) const
Definition: pzmatrix.h:253
TPZSwelling::fLambda
STATE fLambda
Compression modulus [N/mm^2].
Definition: swelling.h:32
TPZVec::NElements
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
rdt.values
def values
Definition: rdt.py:119
TPZSwelling::gState
static TPZFMatrix< STATE > gState
Variables which holds the state variables used in the check convergence procedure.
Definition: swelling.h:279
main
int main()
Definition: benchadolc.cc:42
checkconv.h
Contains the implementation of the CheckConvergence function.
TPZSwelling::ContributeBC
virtual void ContributeBC(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override
It computes a contribution to the stiffness matrix and load vector at one BC integration point...
Definition: swelling.cpp:128
swelling.h
Contains the TPZSwelling class which implements a numerical model of swelling material coupling flow...
TPZMaterialData::sol
TPZSolVec sol
vector of the solutions at the integration point
Definition: pzmaterialdata.h:77
TPZSwelling::gVMinus
static STATE gVMinus
Molar volume anions [mm^3/mmol].
Definition: swelling.h:85
TPZSwelling::NSolutionVariables
virtual int NSolutionVariables(int var) override
Returns the number of solution variables associated with a variable index.
Definition: swelling.cpp:115
TPZMatrix::SolveDirect
virtual int SolveDirect(TPZFMatrix< TVar > &F, const DecomposeType dt, std::list< int64_t > &singular)
Solves the linear system using Direct methods.
Definition: pzmatrix.cpp:710
PZError
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15
TPZBndCond::Material
TPZMaterial * Material() const
Definition: pzbndcond.h:263
TPZSwelling::Residual
void Residual(TPZFMatrix< STATE > &res, int cases)
Computes the residual for the given state variable.
TPZSwelling::NResidual
void NResidual(TPZVec< STATE > &mu, STATE ksi, STATE pressure, TPZVec< STATE > &N, TPZFMatrix< STATE > &res, TPZFMatrix< STATE > &tangent)
Computes the residual and tangent vector of the system of equations which determines N...
Generated on Sun Aug 16 2020 11:17:03 for NeoPZ by   doxygen 1.8.13