Trilinos based (stochastic) FEM solvers
manufacturedSolution.hpp
Go to the documentation of this file.
1 /*
2 Brian Staber (brian.staber@gmail.com)
3 */
4 
5 #ifndef MANUFACTUREDSOLUTION_HPP
6 #define MANUFACTUREDSOLUTION_HPP
7 
8 #include "tensor_calculus.hpp"
9 #include "newtonRaphson.hpp"
10 #include "distributenrldata.hpp"
11 #include "compressibleHyperelasticity.hpp"
12 
13 class manufacturedSolution : public compressibleHyperelasticity
14 {
15 public:
16 
17  double mu1, mu2, mu3, mu4, mu5, mu;
18  double trm, beta3, beta4, beta5;
19  double ptrmbeta4, ptrmbeta5;
20  double plyagl, cos_plyagl, sin_plyagl, topcoord;
21 
22  manufacturedSolution(Epetra_Comm & comm, Teuchos::ParameterList & Parameters, std::string & mesh_file){
23  Mesh = new mesh(comm, mesh_file, 1.0);
24  Comm = Mesh->Comm;
25 
26  StandardMap = new Epetra_Map(-1,3*Mesh->n_local_nodes_without_ghosts,&Mesh->local_dof_without_ghosts[0],0,*Comm);
27  OverlapMap = new Epetra_Map(-1,3*Mesh->n_local_nodes,&Mesh->local_dof[0],0,*Comm);
28  ImportToOverlapMap = new Epetra_Import(*OverlapMap,*StandardMap);
29  create_FECrsGraph();
30  topcoord = 25.0;
31  setup_dirichlet_conditions();
32  }
33 
34  ~manufacturedSolution(){
35  }
36 
37  void set_parameters(Epetra_SerialDenseVector & x, double & angle){
38  mu1 = x(0); mu2 = x(1); mu3 = x(2); mu4 = x(3); mu5 = x(4);
39  beta3 = -0.5; beta4 = x(5); beta5 = x(6);
40  mu = 2.0*mu1 + 4.0*mu2 + 2.0*mu3;
41  trm = mu4 + 2.0*mu5;
42  ptrmbeta4 = std::pow(trm,beta4);
43  ptrmbeta5 = std::pow(trm,beta5);
44  plyagl = angle;
45  cos_plyagl = std::cos(plyagl);
46  sin_plyagl = std::sin(plyagl);
47  }
48 
49  void get_matrix_and_rhs(Epetra_Vector & x, Epetra_FECrsMatrix & K, Epetra_FEVector & F){
50  assembleMixedDirichletNeumann_inhomogeneousForcing(x,K,F);
51  }
52 
53  void setup_dirichlet_conditions(){
54  n_bc_dof = 0;
55  double y;
56  unsigned int node;
57  for (unsigned int i=0; i<Mesh->n_local_nodes_without_ghosts; ++i){
58  node = Mesh->local_nodes[i];
59  y = Mesh->nodes_coord[3*node+1];
60  if(y==0.0){
61  n_bc_dof+=3;
62  }
63  if(y==topcoord){
64  n_bc_dof+=3;
65  }
66  }
67  int indbc = 0;
68  dof_on_boundary = new int [n_bc_dof];
69  for (unsigned int inode=0; inode<Mesh->n_local_nodes_without_ghosts; ++inode){
70  node = Mesh->local_nodes[inode];
71  y = Mesh->nodes_coord[3*node+1];
72  if (y==0.0){
73  dof_on_boundary[indbc+0] = 3*inode+0;
74  dof_on_boundary[indbc+1] = 3*inode+1;
75  dof_on_boundary[indbc+2] = 3*inode+2;
76  indbc+=3;
77  }
78  if (y==topcoord){
79  dof_on_boundary[indbc+0] = 3*inode+0;
80  dof_on_boundary[indbc+1] = 3*inode+1;
81  dof_on_boundary[indbc+2] = 3*inode+2;
82  indbc+=3;
83  }
84  }
85  }
86 
87  void apply_dirichlet_conditions(Epetra_FECrsMatrix & K, Epetra_FEVector & F, double & displacement){
88  Epetra_MultiVector v(*StandardMap,true);
89  Epetra_SerialDenseVector u(3);
90  v.PutScalar(0.0);
91  double x,y,z;
92  int node;
93  for (unsigned int inode=0; inode<Mesh->n_local_nodes_without_ghosts; ++inode){
94  node = Mesh->local_nodes[inode];
95  x = Mesh->nodes_coord[3*node+0];
96  y = Mesh->nodes_coord[3*node+1];
97  z = Mesh->nodes_coord[3*node+2];
98  if (y==0.0){
99  u = getManufacturedSolution(x,y,z);
100  u.Scale(displacement);
101  v[0][StandardMap->LID(3*node+0)] = u(0);
102  v[0][StandardMap->LID(3*node+1)] = u(1);
103  v[0][StandardMap->LID(3*node+2)] = u(2);
104  }
105  if (y==topcoord){
106  u = getManufacturedSolution(x,y,z);
107  u.Scale(displacement);
108  v[0][StandardMap->LID(3*node+0)] = u(0);
109  v[0][StandardMap->LID(3*node+1)] = u(1);
110  v[0][StandardMap->LID(3*node+2)] = u(2);
111  }
112  }
113  Epetra_MultiVector rhs_dir(*StandardMap,true);
114  K.Apply(v,rhs_dir);
115  F.Update(-1.0,rhs_dir,1.0);
116  for (unsigned int inode=0; inode<Mesh->n_local_nodes_without_ghosts; ++inode){
117  node = Mesh->local_nodes[inode];
118  y = Mesh->nodes_coord[3*node+1];
119  if (y==0.0){
120  F[0][StandardMap->LID(3*node+0)] = v[0][StandardMap->LID(3*node+0)];
121  F[0][StandardMap->LID(3*node+1)] = v[0][StandardMap->LID(3*node+1)];
122  F[0][StandardMap->LID(3*node+2)] = v[0][StandardMap->LID(3*node+2)];
123  }
124  if (y==topcoord){
125  F[0][StandardMap->LID(3*node+0)] = v[0][StandardMap->LID(3*node+0)];
126  F[0][StandardMap->LID(3*node+1)] = v[0][StandardMap->LID(3*node+1)];
127  F[0][StandardMap->LID(3*node+2)] = v[0][StandardMap->LID(3*node+2)];
128  }
129  }
130  ML_Epetra::Apply_OAZToMatrix(dof_on_boundary,n_bc_dof,K);
131  }
132 
133  Epetra_SerialDenseVector get_neumannBc(Epetra_SerialDenseMatrix & matrix_X, Epetra_SerialDenseMatrix & xg, unsigned int & gp){
134  Epetra_SerialDenseVector h(3), normal(3);
135  Epetra_SerialDenseMatrix piola(3,3), d_shape_functions(Mesh->face_type,2), dxi_matrix_x(3,2);
136  piola = getManufacturedPiola(xg(0,gp),xg(1,gp),xg(2,gp));
137  for (unsigned int inode=0; inode<Mesh->face_type; ++inode){
138  d_shape_functions(inode,0) = Mesh->D1_N_faces(gp,inode);
139  d_shape_functions(inode,1) = Mesh->D2_N_faces(gp,inode);
140  }
141  dxi_matrix_x.Multiply('N','N',1.0,matrix_X,d_shape_functions,0.0);
142  normal(0) = dxi_matrix_x(1,0)*dxi_matrix_x(2,1) - dxi_matrix_x(2,0)*dxi_matrix_x(1,1);
143  normal(1) = dxi_matrix_x(2,0)*dxi_matrix_x(0,1) - dxi_matrix_x(0,0)*dxi_matrix_x(2,1);
144  normal(2) = dxi_matrix_x(0,0)*dxi_matrix_x(1,1) - dxi_matrix_x(1,0)*dxi_matrix_x(0,1);
145  if (xg(2,gp)==0.0){
146  normal.Scale(-1.0);
147  }
148  h.Multiply('N','N',1.0,piola,normal,0.0);
149  return h;
150  }
151 
152  Epetra_SerialDenseVector get_forcing(double & x1, double & x2, double & x3, unsigned int & e_lid, unsigned int & gp){
153  Epetra_SerialDenseVector f(3);
154  f = manufacturedForcing(x1,x2,x3);
155  return f;
156  }
157 
158  void get_material_parameters(unsigned int & e_lid, unsigned int & gp){
159  //std::cout << "**Err: Not using that method in this example!\n";
160  }
161 
162  void get_constitutive_tensors(Epetra_SerialDenseMatrix & deformation_gradient, Epetra_SerialDenseVector & piola_stress, Epetra_SerialDenseMatrix & tangent_piola){
163 
164  double det = deformation_gradient(0,0)*deformation_gradient(1,1)*deformation_gradient(2,2)-deformation_gradient(0,0)*deformation_gradient(1,2)*deformation_gradient(2,1)-deformation_gradient(0,1)*deformation_gradient(1,0)*deformation_gradient(2,2)+deformation_gradient(0,1)*deformation_gradient(1,2)*deformation_gradient(2,0)+deformation_gradient(0,2)*deformation_gradient(1,0)*deformation_gradient(2,1)-deformation_gradient(0,2)*deformation_gradient(1,1)*deformation_gradient(2,0);
165 
166  Epetra_SerialDenseMatrix C(3,3);
167  Epetra_SerialDenseMatrix CC(3,3);
168  Epetra_SerialDenseMatrix ddJ5(6,6);
169  Epetra_SerialDenseVector LML(6);
170  Epetra_SerialDenseVector eye(6);
171  Epetra_SerialDenseVector dJ5(6);
172  Epetra_SerialDenseVector M(6);
173  Epetra_SerialDenseVector L(6);
174  Epetra_SerialDenseVector c(6);
175 
176  C.Multiply('T','N',1.0,deformation_gradient,deformation_gradient,0.0);
177  CC.Multiply('N','N',1.0,C,C,0.0);
178 
179  c(0) = C(0,0); c(1) = C(1,1); c(2) = C(2,2); c(3) = C(1,2); c(4) = C(0,2); c(5) = C(0,1);
180 
181  L(0) = (1.0/(det*det))*(C(1,1)*C(2,2)-C(1,2)*C(2,1));
182  L(1) = (1.0/(det*det))*(C(0,0)*C(2,2)-C(0,2)*C(2,0));
183  L(2) = (1.0/(det*det))*(C(0,0)*C(1,1)-C(0,1)*C(1,0));
184  L(3) = (1.0/(det*det))*(C(0,2)*C(1,0)-C(0,0)*C(1,2));
185  L(4) = (1.0/(det*det))*(C(0,1)*C(1,2)-C(0,2)*C(1,1));
186  L(5) = (1.0/(det*det))*(C(0,2)*C(2,1)-C(0,1)*C(2,2));
187 
188  M(0) = mu4*sin_plyagl*sin_plyagl+mu5*cos_plyagl*cos_plyagl;
189  M(1) = mu4*cos_plyagl*cos_plyagl+mu5*sin_plyagl*sin_plyagl;
190  M(2) = mu5;
191  M(3) = 0.0;
192  M(4) = 0.0;
193  M(5) = (mu4-mu5)*cos_plyagl*sin_plyagl;
194 
195  LML(0) = M(0)*L(0)*L(0)+2.0*M(4)*L(0)*L(4)+2.0*M(5)*L(0)*L(5)+M(2)*L(4)*L(4)+2.0*M(3)*L(4)*L(5)+M(1)*L(5)*L(5);
196  LML(1) = M(1)*L(1)*L(1)+2.0*M(4)*L(1)*L(3)+2.0*M(5)*L(1)*L(5)+M(2)*L(3)*L(3)+2.0*M(4)*L(3)*L(5)+M(0)*L(5)*L(5);
197  LML(2) = M(2)*L(2)*L(2)+2.0*M(3)*L(2)*L(3)+2.0*M(4)*L(2)*L(4)+M(1)*L(3)*L(3)+2.0*M(5)*L(3)*L(4)+M(0)*L(4)*L(4);
198  LML(3) = L(2)*(L(1)*M(3)+L(3)*M(2)+L(5)*M(4))+L(3)*(L(1)*M(1)+L(3)*M(3)+L(5)*M(5))+L(4)*(L(5)*M(0)+L(1)*M(5)+L(3)*M(4));
199  LML(4) = L(2)*(L(0)*M(4)+L(4)*M(2)+L(5)*M(3))+L(3)*(L(0)*M(5)+L(5)*M(1)+L(4)*M(3))+L(4)*(L(0)*M(0)+L(4)*M(4)+L(5)*M(5));
200  LML(5) = L(1)*(L(0)*M(5)+L(5)*M(1)+L(4)*M(3))+L(3)*(L(0)*M(4)+L(4)*M(2)+L(5)*M(3))+L(5)*(L(0)*M(0)+L(4)*M(4)+L(5)*M(5));
201 
202  eye(0) = 1.0; eye(1) = 1.0; eye(2) = 1.0; eye(3) = 0.0; eye(4) = 0.0; eye(5) = 0.0;
203 
204  double I1 = C(0,0) + C(1,1) + C(2,2);
205  double II1 = C(0,0)*C(0,0) + C(1,1)*C(1,1) + C(2,2)*C(2,2) + 2.0*C(1,2)*C(1,2) + 2.0*C(0,2)*C(0,2) + 2.0*C(0,1)*C(0,1);
206  double I2 = (1.0/2.0)*(I1*I1-II1);
207  double I3 = det*det;
208  double I4 = C(0,0)*M(0) + C(1,1)*M(1) + C(2,2)*M(2) + 2.0*C(0,1)*M(5) + 2.0*C(0,2)*M(4) + 2.0*C(1,2)*M(3);
209  double I5 = CC(0,0)*M(0) + CC(1,1)*M(1) + CC(2,2)*M(2) + 2.0*CC(0,1)*M(5) + 2.0*CC(0,2)*M(4) + 2.0*CC(1,2)*M(3);
210  double J5 = I5 - I1*I4 + I2*trm;
211  double pI3 = std::pow(I3,-beta3);
212  double pI4 = std::pow(I4,beta4);
213  double pJ5 = std::pow(J5,beta5);
214 
215  for (unsigned int i=0; i<6; ++i){
216  dJ5(i) = J5*L(i) - I3*LML(i);
217  piola_stress(i) = 2.0*mu1*eye(i) + 2.0*mu2*(I1*eye(i)-c(i)) + (2.0*mu3*det*det-mu)*L(i) + (2.0/ptrmbeta4)*pI4*M(i) + (2.0/ptrmbeta5)*pJ5*dJ5(i) - 2.0*trm*pI3*L(i);
218  }
219 
220  double scalarAB = 4.0*mu2;
221  tensor_product(scalarAB,eye,eye,tangent_piola,0.0);
222 
223  scalarAB = -4.0*mu2;
224  sym_tensor_product(scalarAB,eye,eye,tangent_piola,1.0);
225 
226  scalarAB = 4.0*mu3*det*det;
227  tensor_product(scalarAB,L,L,tangent_piola,1.0);
228 
229  scalarAB = -4.0*mu3*det*det+2.0*mu;
230  sym_tensor_product(scalarAB,L,L,tangent_piola,1.0);
231 
232  scalarAB = J5;
233  tensor_product(scalarAB,L,L,ddJ5,0.0);
234  scalarAB = -J5;
235  sym_tensor_product(scalarAB,L,L,ddJ5,1.0);
236  scalarAB = -I3;
237  tensor_product(scalarAB,L,LML,ddJ5,1.0);
238  tensor_product(scalarAB,LML,L,ddJ5,1.0);
239  scalarAB = I3;
240  sym_tensor_product(scalarAB,L,LML,ddJ5,1.0);
241  sym_tensor_product(scalarAB,LML,L,ddJ5,1.0);
242 
243  scalarAB = 4.0*beta4*pI4/(I4*ptrmbeta4);
244  tensor_product(scalarAB,M,M,tangent_piola,1.0);
245 
246  scalarAB = 4.0*beta5*pJ5/(J5*ptrmbeta5);
247  tensor_product(scalarAB,dJ5,dJ5,tangent_piola,1.0);
248 
249  scalarAB = 4.0*trm*beta3*pI3;
250  tensor_product(scalarAB,L,L,tangent_piola,1.0);
251  scalarAB = 4.0*trm*pI3;
252  sym_tensor_product(scalarAB,L,L,tangent_piola,1.0);
253 
254  ddJ5.Scale(4.0*pJ5/ptrmbeta5);
255  tangent_piola += ddJ5;
256  }
257 
258 
259  Epetra_SerialDenseVector getManufacturedSolution(double & x1, double & x2, double & x3){
260  Epetra_SerialDenseVector u(3);
261  double c1 = 2.0e-4; double c2 = 1.0e-4; double c3 = 2.0e-4;
262  u(0) = 0.0;
263  u(1) = x2*x2*1.1/(25.0*25.0);
264  u(2) = 0.0;
265  /*u(0) = -c1*x2*(topcoord-x1)*(topcoord-x2);
266  u(1) = c2*x2*(topcoord-x2);
267  u(2) = std::sin(c1*x3);*/
268  return u;
269  }
270 
271  Epetra_SerialDenseMatrix getManufacturedPiola(double & x1, double & x2, double & x3){
272 
273  Epetra_SerialDenseMatrix F(3,3), C(3,3), CC(3,3), ML(3,3), LML(3,3), L(3,3), M(3,3), eye(3,3), S(3,3), P(3,3);
274  double c1 = 2.0e-4; double c2 = 1.0e-4; double c3 = 2.0e-4;
275 
276  /*F(0,0) = c1*x2*(topcoord-x2)+1.0; F(0,1) = c1*x2*(topcoord-x1)-c1*(topcoord-x1)*(topcoord-x2); F(0,2) = 0.0;
277  F(1,0) = 0.0; F(1,1) = c2*(topcoord-x2)-c2*x2+1.0; F(1,2) = 0.0;
278  F(2,0) = 0.0; F(2,1) = 0.0; F(2,2) = c1*std::cos(c1*x3)+1.0;*/
279  F(0,0) = 1.0; F(0,1) = 0.0; F(0,2) = 0.0;
280  F(1,0) = 0.0; F(1,1) = 1.0+2.0*1.1*x2/(25.0*25.0); F(1,2) = 0.0;
281  F(2,0) = 0.0; F(2,1) = 0.0; F(2,2) = 1.0;
282  double det = F(0,0)*F(1,1)*F(2,2)-F(0,0)*F(1,2)*F(2,1)-F(0,1)*F(1,0)*F(2,2)+F(0,1)*F(1,2)*F(2,0)+F(0,2)*F(1,0)*F(2,1)-F(0,2)*F(1,1)*F(2,0);
283 
284  M(0,0) = mu4*sin_plyagl*sin_plyagl+mu5*cos_plyagl*cos_plyagl;
285  M(1,1) = mu4*cos_plyagl*cos_plyagl+mu5*sin_plyagl*sin_plyagl;
286  M(2,2) = mu5;
287  M(1,2) = 0.0; M(2,1) = 0.0;
288  M(0,2) = 0.0; M(2,0) = 0.0;
289  M(0,1) = (mu4-mu5)*cos_plyagl*sin_plyagl; M(1,0) = M(0,1);
290 
291  C.Multiply('T','N',1.0,F,F,0.0);
292  CC.Multiply('N','N',1.0,C,C,0.0);
293 
294  L(0,0) = (1.0/(det*det))*(C(1,1)*C(2,2)-C(1,2)*C(2,1));
295  L(1,1) = (1.0/(det*det))*(C(0,0)*C(2,2)-C(0,2)*C(2,0));
296  L(2,2) = (1.0/(det*det))*(C(0,0)*C(1,1)-C(0,1)*C(1,0));
297  L(1,2) = (1.0/(det*det))*(C(0,2)*C(1,0)-C(0,0)*C(1,2)); L(2,1) = L(1,2);
298  L(0,2) = (1.0/(det*det))*(C(0,1)*C(1,2)-C(0,2)*C(1,1)); L(2,0) = L(0,2);
299  L(0,1) = (1.0/(det*det))*(C(0,2)*C(2,1)-C(0,1)*C(2,2)); L(1,0) = L(0,1);
300 
301  ML.Multiply('N','N',1.0,M,L,0.0);
302  LML.Multiply('N','N',1.0,L,ML,0.0);
303 
304  double I1 = C(0,0) + C(1,1) + C(2,2);
305  double II1 = C(0,0)*C(0,0) + C(1,1)*C(1,1) + C(2,2)*C(2,2) + 2.0*C(1,2)*C(1,2) + 2.0*C(0,2)*C(0,2) + 2.0*C(0,1)*C(0,1);
306  double I2 = (1.0/2.0)*(I1*I1-II1);
307  double I3 = det*det;
308  double I4 = C(0,0)*M(0,0) + C(1,1)*M(1,1) + C(2,2)*M(2,2) + 2.0*C(0,1)*M(0,1) + 2.0*C(0,2)*M(0,2) + 2.0*C(1,2)*M(1,2);
309  double I5 = CC(0,0)*M(0,0) + CC(1,1)*M(1,1) + CC(2,2)*M(2,2) + 2.0*CC(0,1)*M(0,1) + 2.0*CC(0,2)*M(0,2) + 2.0*CC(1,2)*M(1,2);
310  double J5 = I5 - I1*I4 + I2*trm;
311  double pI3 = std::pow(I3,-beta3);
312  double pI4 = std::pow(I4,beta4);
313  double pJ5 = std::pow(J5,beta5);
314 
315  eye(0,0) = 1.0; eye(0,1) = 0.0; eye(0,2) = 0.0;
316  eye(1,0) = 0.0; eye(1,1) = 1.0; eye(1,2) = 0.0;
317  eye(2,0) = 0.0; eye(2,1) = 0.0; eye(2,2) = 1.0;
318 
319  for (unsigned int i=0; i<3; ++i){
320  for (unsigned int j=0; j<3; ++j){
321  S(i,j) = 2.0*mu1*eye(i,j) + 2.0*mu2*(I1*eye(i,j)-C(i,j)) + (2.0*mu3*det*det-mu)*L(i,j) + (2.0/ptrmbeta4)*pI4*M(i,j) + (2.0/ptrmbeta5)*pJ5*(J5*L(i,j)-I3*LML(i,j)) - 2.0*trm*pI3*L(i,j);
322  }
323  }
324  P.Multiply('N','N',1.0,F,S,0.0);
325  return P;
326  }
327 
328  Epetra_SerialDenseVector manufacturedForcing(double & x1, double & x2, double & x3){
329  Epetra_SerialDenseVector f(3);
330  /*Epetra_SerialDenseVector xf(3), xb(3);
331  Epetra_SerialDenseMatrix Pf(3,3), Pb(3,3);
332  double h = 1.0e-5;
333  xf(0) = x1+h; xf(1) = x2; xf(2) = x3;
334  xb(0) = x1-h; xb(1) = x2; xb(2) = x3;
335  Pf = getManufacturedPiola(xf(0),xf(1),xf(2));
336  Pb = getManufacturedPiola(xb(0),xb(1),xb(2));
337  f(0) = (Pf(0,0)-Pb(0,0))/(2.0*h);
338  f(1) = (Pf(1,0)-Pb(1,0))/(2.0*h);
339  f(2) = (Pf(2,0)-Pb(2,0))/(2.0*h);
340 
341  xf(0) = x1; xf(1) = x2+h; xf(2) = x3;
342  xb(0) = x1; xb(1) = x2-h; xb(2) = x3;
343  Pf = getManufacturedPiola(xf(0),xf(1),xf(2));
344  Pb = getManufacturedPiola(xb(0),xb(1),xb(2));
345  f(0) += (Pf(0,1)-Pb(0,1))/(2.0*h);
346  f(1) += (Pf(1,1)-Pb(1,1))/(2.0*h);
347  f(2) += (Pf(2,1)-Pb(2,1))/(2.0*h);
348 
349  xf(0) = x1; xf(1) = x2; xf(2) = x3+h;
350  xb(0) = x1; xb(1) = x2; xb(2) = x3-h;
351  Pf = getManufacturedPiola(xf(0),xf(1),xf(2));
352  Pb = getManufacturedPiola(xb(0),xb(1),xb(2));
353  f(0) += (Pf(0,2)-Pb(0,2))/(2.0*h);
354  f(1) += (Pf(1,2)-Pb(1,2))/(2.0*h);
355  f(2) += (Pf(2,2)-Pb(2,2))/(2.0*h);*/
356 
357  double M11 = mu4*sin_plyagl*sin_plyagl+mu5*cos_plyagl*cos_plyagl;
358  double M22 = mu4*cos_plyagl*cos_plyagl+mu5*sin_plyagl*sin_plyagl;
359  double M33 = mu5;
360  double M23 = 0.0; double M32 = 0.0;
361  double M13 = 0.0; double M31 = 0.0;
362  double M12 = (mu4-mu5)*cos_plyagl*sin_plyagl; double M21 = M12;
363 
364  double ud = 1.1/(25.0*25.0);
365 
366  f(0) = (8*M12*M22*beta4*ud*(2*ud*x2+1)*std::pow(4*M22*ud*ud*x2*x2+4*M22*ud*x2+M11+M22+M33,beta4-1))/std::pow(M11+M22+M33,beta4) - (8*M12*beta5*ud*(2*ud*x2+1)*(M11+M33)*std::pow(M11+M22+M33 + 4*M11*ud*ud*x2*x2 + 4*M33*ud*ud*x2*x2 + 4*M11*ud*x2 + 4*M33*ud*x2,beta5-1))/std::pow(M11+M22+M33,beta5);
367  f(1) = 2*ud*(2*mu1 + 4*mu2 - (- 8*mu3*ud*ud*x2*x2 - 8*mu3*ud*x2 + 2*mu1 + 4*mu2)/(2*ud*x2 + 1)*(2*ud*x2 + 1) + (2*M22*std::pow(4*M22*ud*ud*x2*x2 + 4*M22*ud*x2 + M11 + M22 + M33,beta4))/std::pow(M11 + M22 + M33,beta4) + (2*(M11 + M33)*std::pow(M11 + M22 + M33 + 4*M11*ud*ud*x2*x2 + 4*M33*ud*ud*x2*x2 + 4*M11*ud*x2 + 4*M33*ud*x2,beta5))/std::pow(M11 + M22 + M33,beta5) - (sign(2*ud*x2 + 1)*(2*M11 + 2*M22 + 2*M33))/(2*ud*x2 + 1)) + (2*ud*x2 + 1)*((4*ud*(- 8*mu3*ud*ud*x2*x2 - 8*mu3*ud*x2 + 2*mu1 + 4*mu2))/std::pow(2*ud*x2 + 1,3) + (8*mu3*ud)/(2*ud*x2 + 1) + (2*ud*sign(2*ud*x2 + 1)*(2*M11 + 2*M22 + 2*M33))/(2*ud*x2 + 1)*(2*ud*x2 + 1) + (8*M22*M22*beta4*ud*(2*ud*x2 + 1)*std::pow(4*M22*ud*ud*x2*x2 + 4*M22*ud*x2 + M11 + M22 + M33,beta4 - 1))/std::pow(M11 + M22 + M33,beta4) + (8*beta5*ud*(2*ud*x2 + 1)*(M11 + M33)*(M11 + M33)*std::pow(M11 + M22 + M33 + 4*M11*ud*ud*x2*x2 + 4*M33*ud*ud*x2*x2 + 4*M11*ud*x2 + 4*M33*ud*x2,beta5 - 1))/std::pow(M11 + M22 + M33,beta5));
368  f(2) = 0.0;
369 
370  f.Scale(-1.0);
371  return f;
372  }
373 
374  template <typename T>
375  int sign (const T &val) { return (val > 0) - (val < 0); }
376 
377  void get_material_parameters_for_recover(unsigned int & e_lid){
378  std::cout << "**Err: Not using that method in this example!\n";
379  }
380 
381  void get_stress_for_recover(Epetra_SerialDenseMatrix & deformation_gradient, double & det, Epetra_SerialDenseMatrix & piola_stress){
382  std::cout << "**Err: Not using that method in this example!\n";
383  }
384 
385  double errorL2(Epetra_Vector & uStandardMap){
386 
387  Epetra_Vector u(*OverlapMap);
388  u.Import(uStandardMap, *ImportToOverlapMap, Insert);
389  double totalError = 0.0;
390  double error = 0.0;
391  double normVH;
392  double gauss_weight;
393  int n_gauss_points = Mesh->n_gauss_cells;
394  int e_gid, node;
395 
396  //Epetra_SerialDenseVector epsilon(6);
397  //Epetra_SerialDenseMatrix matrix_B(6,3*Mesh->el_type), dx_shape_functions(Mesh->el_type,3);
398  Epetra_SerialDenseVector uExact(3), vH(3);
399  Epetra_SerialDenseMatrix X_I(3,Mesh->el_type), u_I(3,Mesh->el_type);
400  Epetra_SerialDenseMatrix u_G(3,n_gauss_points), x_G(3,n_gauss_points);
401 
402  for (unsigned int e_lid=0; e_lid<Mesh->n_local_cells; ++e_lid){
403  e_gid = Mesh->local_cells[e_lid];
404  for (unsigned int inode=0; inode<Mesh->el_type; ++inode){
405  node = Mesh->cells_nodes[Mesh->el_type*e_gid+inode];
406  X_I(0,inode) = Mesh->nodes_coord[3*node+0];
407  X_I(1,inode) = Mesh->nodes_coord[3*node+1];
408  X_I(2,inode) = Mesh->nodes_coord[3*node+2];
409  u_I(0,inode) = u[OverlapMap->LID(3*node+0)];
410  u_I(1,inode) = u[OverlapMap->LID(3*node+1)];
411  u_I(2,inode) = u[OverlapMap->LID(3*node+2)];
412  }
413  x_G.Multiply('N','N',1.0,X_I,Mesh->N_cells,0.0);
414  u_G.Multiply('N','N',1.0,u_I,Mesh->N_cells,0.0);
415  for (unsigned int gp=0; gp<n_gauss_points; ++gp){
416  gauss_weight = Mesh->gauss_weight_cells(gp);
417  uExact = getManufacturedSolution(x_G(0,gp),x_G(1,gp),x_G(2,gp));
418  vH(0) = uExact(0) - u_G(0,gp);
419  vH(1) = uExact(1) - u_G(1,gp);
420  vH(2) = uExact(2) - u_G(2,gp);
421  normVH = vH.Norm2();
422  error += gauss_weight*normVH*normVH*Mesh->detJac_cells(e_lid,gp);
423  /*for (unsigned int inode=0; inode<Mesh->el_type; ++inode){
424  dx_shape_functions(inode,0) = Mesh->DX_N_cells(gp+n_gauss_points*inode,e_lid);
425  dx_shape_functions(inode,1) = Mesh->DY_N_cells(gp+n_gauss_points*inode,e_lid);
426  dx_shape_functions(inode,2) = Mesh->DZ_N_cells(gp+n_gauss_points*inode,e_lid);
427  }
428  compute_B_matrices(dx_shape_functions,matrix_B);
429  epsilon.Multiply('N','N',1.0,matrix_B,vector_u,0.0);*/
430  }
431  }
432  Comm->SumAll(&error,&totalError,1);
433  totalError = std::sqrt(totalError);
434  return totalError;
435  }
436 
437 };
438 
439 #endif
tensor_calculus.hpp
manufacturedSolution::get_forcing
Epetra_SerialDenseVector get_forcing(double &x1, double &x2, double &x3, unsigned int &e_lid, unsigned int &gp)
Definition: manufacturedSolution.hpp:152
baseClassFEM::Comm
Epetra_Comm * Comm
Definition: baseClassFEM.hpp:15
j
for j
Definition: costFunction.m:42
u
u
Definition: run.m:9
mesh::local_nodes
std::vector< int > local_nodes
Definition: meshpp.hpp:80
tensor_product
void tensor_product(double scalarAB, Epetra_SerialDenseVector &A, Epetra_SerialDenseVector &B, Epetra_SerialDenseMatrix &AoB, double scalarThis)
Definition: tensor_calculus.hpp:10
manufacturedSolution::errorL2
double errorL2(Epetra_Vector &uStandardMap)
Definition: manufacturedSolution.hpp:385
mesh::nodes_coord
std::vector< double > nodes_coord
Definition: meshpp.hpp:77
mesh::local_dof
std::vector< int > local_dof
Definition: meshpp.hpp:80
manufacturedSolution::manufacturedForcing
Epetra_SerialDenseVector manufacturedForcing(double &x1, double &x2, double &x3)
Definition: manufacturedSolution.hpp:328
hyperelasticity::n_bc_dof
unsigned int n_bc_dof
Definition: hyperelasticity.hpp:23
mesh::D2_N_faces
Epetra_SerialDenseMatrix D2_N_faces
Definition: meshpp.hpp:71
mesh::N_cells
Epetra_SerialDenseVector N_cells
Definition: meshpp.hpp:73
mesh::Comm
Epetra_Comm * Comm
Definition: meshpp.hpp:69
mesh::n_local_cells
int n_local_cells
Definition: meshpp.hpp:94
compressibleHyperelasticity::assembleMixedDirichletNeumann_inhomogeneousForcing
void assembleMixedDirichletNeumann_inhomogeneousForcing(Epetra_Vector &x, Epetra_FECrsMatrix &K, Epetra_FEVector &F)
Definition: compressibleHyperelasticity.cpp:59
mesh::face_type
int face_type
Definition: meshpp.hpp:91
mesh::detJac_cells
Epetra_SerialDenseVector detJac_cells
Definition: meshpp.hpp:73
mesh::local_cells
std::vector< int > local_cells
Definition: meshpp.hpp:81
mesh
clc clear all close all mesh
Definition: run.m:5
mesh::n_local_nodes_without_ghosts
int n_local_nodes_without_ghosts
Definition: meshpp.hpp:92
mesh::n_local_nodes
int n_local_nodes
Definition: meshpp.hpp:93
mesh::el_type
int el_type
Definition: meshpp.hpp:90
manufacturedSolution::apply_dirichlet_conditions
void apply_dirichlet_conditions(Epetra_FECrsMatrix &K, Epetra_FEVector &F, double &displacement)
Definition: manufacturedSolution.hpp:87
mesh::cells_nodes
std::vector< int > cells_nodes
Definition: meshpp.hpp:78
L
L
Definition: costFunction.m:66
baseClassFEM::ImportToOverlapMap
Epetra_Import * ImportToOverlapMap
Definition: baseClassFEM.hpp:19
baseClassFEM::StandardMap
Epetra_Map * StandardMap
Definition: baseClassFEM.hpp:18
mesh::D1_N_faces
Epetra_SerialDenseMatrix D1_N_faces
Definition: meshpp.hpp:71
manufacturedSolution::getManufacturedPiola
Epetra_SerialDenseMatrix getManufacturedPiola(double &x1, double &x2, double &x3)
Definition: manufacturedSolution.hpp:271
mesh::n_gauss_cells
unsigned int n_gauss_cells
Definition: meshpp.hpp:98
manufacturedSolution::get_matrix_and_rhs
void get_matrix_and_rhs(Epetra_Vector &x, Epetra_FECrsMatrix &K, Epetra_FEVector &F)
Definition: manufacturedSolution.hpp:49
sym_tensor_product
void sym_tensor_product(double scalarAB, Epetra_SerialDenseVector &A, Epetra_SerialDenseVector &B, Epetra_SerialDenseMatrix &AoB, double scalarThis)
Definition: tensor_calculus.hpp:16
manufacturedSolution::set_parameters
void set_parameters(Epetra_SerialDenseVector &x, double &angle)
Definition: manufacturedSolution.hpp:37
manufacturedSolution::get_stress_for_recover
void get_stress_for_recover(Epetra_SerialDenseMatrix &deformation_gradient, double &det, Epetra_SerialDenseMatrix &piola_stress)
Definition: manufacturedSolution.hpp:381
mesh::local_dof_without_ghosts
std::vector< int > local_dof_without_ghosts
Definition: meshpp.hpp:79
manufacturedSolution::get_material_parameters_for_recover
void get_material_parameters_for_recover(unsigned int &e_lid)
Definition: manufacturedSolution.hpp:377
manufacturedSolution::get_material_parameters
void get_material_parameters(unsigned int &e_lid, unsigned int &gp)
Definition: manufacturedSolution.hpp:158
tri3::d_shape_functions
void d_shape_functions(Epetra_SerialDenseMatrix &D, double &xi, double &eta)
Definition: fepp.cpp:12
i
for i
Definition: costFunction.m:38
baseClassFEM::OverlapMap
Epetra_Map * OverlapMap
Definition: baseClassFEM.hpp:17
manufacturedSolution::get_constitutive_tensors
void get_constitutive_tensors(Epetra_SerialDenseMatrix &deformation_gradient, Epetra_SerialDenseVector &piola_stress, Epetra_SerialDenseMatrix &tangent_piola)
Definition: manufacturedSolution.hpp:162
manufacturedSolution::getManufacturedSolution
Epetra_SerialDenseVector getManufacturedSolution(double &x1, double &x2, double &x3)
Definition: manufacturedSolution.hpp:259
manufacturedSolution::get_neumannBc
Epetra_SerialDenseVector get_neumannBc(Epetra_SerialDenseMatrix &matrix_X, Epetra_SerialDenseMatrix &xg, unsigned int &gp)
Definition: manufacturedSolution.hpp:133
manufacturedSolution::manufacturedSolution
manufacturedSolution(Epetra_Comm &comm, Teuchos::ParameterList &Parameters, std::string &mesh_file)
Definition: manufacturedSolution.hpp:22
mesh::gauss_weight_cells
Epetra_SerialDenseVector gauss_weight_cells
Definition: meshpp.hpp:99