Trilinos based (stochastic) FEM solvers
compressibleHyperelasticity_linearPatchTest.hpp
Go to the documentation of this file.
1 /*
2 Brian Staber (brian.staber@gmail.com)
3 */
4 
5 #ifndef COMPRESSIBLEHYPERELASTICITY_LINEARPATCHTEST_HPP
6 #define COMPRESSIBLEHYPERELASTICITY_LINEARPATCHTEST_HPP
7 
8 #include "tensor_calculus.hpp"
10 #include "newtonRaphson.hpp"
11 
13 {
14 public:
15 
16  double mu1, mu2, mu3, mu4, mu5, mu;
17  double trm, beta3, beta4, beta5;
20 
21  compressibleHyperelasticity_linearPatchTest(Epetra_Comm & comm, Teuchos::ParameterList & Parameters){
22 
23  std::string mesh_file = Teuchos::getParameter<std::string>(Parameters.sublist("Mesh"), "mesh_file");
24  Mesh = new mesh(comm, mesh_file, 1000.0);
25  Comm = Mesh->Comm;
26 
27  StandardMap = new Epetra_Map(-1,3*Mesh->n_local_nodes_without_ghosts,&Mesh->local_dof_without_ghosts[0],0,*Comm);
28  OverlapMap = new Epetra_Map(-1,3*Mesh->n_local_nodes,&Mesh->local_dof[0],0,*Comm);
29  ImportToOverlapMap = new Epetra_Import(*OverlapMap,*StandardMap);
30 
33  }
34 
36  }
37 
38  void set_parameters(Epetra_SerialDenseVector & x, double & angle){
39  mu1 = x(0);
40  mu2 = x(1);
41  mu3 = x(2);
42  mu4 = x(3);
43  mu5 = x(4);
44  beta3 = -0.5;
45  beta4 = x(5);
46  beta5 = x(6);
47  mu = 2.0*mu1 + 4.0*mu2 + 2.0*mu3;
48  trm = mu4 + 2.0*mu5;
49  ptrmbeta4 = std::pow(trm,beta4);
50  ptrmbeta5 = std::pow(trm,beta5);
51  plyagl = angle;
52  cos_plyagl = std::cos(plyagl);
53  sin_plyagl = std::sin(plyagl);
54  }
55 
56  Epetra_SerialDenseVector exactSolution(double & x1, double & x2, double & x3){
57  Epetra_SerialDenseVector u(3);
58  u(0) = 0.2*x1;
59  u(1) = 0.2*x2;
60  u(2) = x3;
61  return u;
62  }
63 
64  double errorL2(Epetra_Vector & uStandardMap){
65 
66  Epetra_Vector u(*OverlapMap);
67  u.Import(uStandardMap, *ImportToOverlapMap, Insert);
68  double totalError;
69  double error = 0.0;
70  double normVH;
71  double gauss_weight;
72  int n_gauss_points = Mesh->n_gauss_cells;
73  int e_gid, node;
74 
75  Epetra_SerialDenseVector uExact(3), vH(3);
76  Epetra_SerialDenseMatrix dx_shape_functions(Mesh->el_type,3), X_I(3,Mesh->el_type), u_I(3,Mesh->el_type);
77  Epetra_SerialDenseMatrix u_G(3,n_gauss_points), x_G(3,n_gauss_points);
78 
79  for (unsigned int e_lid=0; e_lid<Mesh->n_local_cells; ++e_lid){
80  e_gid = Mesh->local_cells[e_lid];
81  for (unsigned int inode=0; inode<Mesh->el_type; ++inode){
82  node = Mesh->cells_nodes[Mesh->el_type*e_gid+inode];
83  X_I(0,inode) = Mesh->nodes_coord[3*node+0];
84  X_I(1,inode) = Mesh->nodes_coord[3*node+1];
85  X_I(2,inode) = Mesh->nodes_coord[3*node+2];
86  u_I(0,inode) = u[OverlapMap->LID(3*node+0)];
87  u_I(1,inode) = u[OverlapMap->LID(3*node+1)];
88  u_I(2,inode) = u[OverlapMap->LID(3*node+2)];
89  }
90  x_G.Multiply('N','N',1.0,X_I,Mesh->N_cells,0.0);
91  u_G.Multiply('N','N',1.0,u_I,Mesh->N_cells,0.0);
92  for (unsigned int gp=0; gp<n_gauss_points; ++gp){
93  gauss_weight = Mesh->gauss_weight_cells(gp);
94  uExact = exactSolution(x_G(0,gp),x_G(1,gp),x_G(2,gp));
95  vH(0) = uExact(0) - u_G(0,gp);
96  vH(1) = uExact(1) - u_G(1,gp);
97  vH(2) = uExact(2) - u_G(2,gp);
98  normVH = vH.Norm2();
99  error += gauss_weight*normVH*normVH*Mesh->detJac_cells(e_lid,gp);
100  }
101  }
102  Comm->SumAll(&error,&totalError,1);
103  totalError = std::sqrt(totalError);
104  return totalError;
105  }
106 
107  void get_matrix_and_rhs(Epetra_Vector & x, Epetra_FECrsMatrix & K, Epetra_FEVector & F){
109  }
110 
112  n_bc_dof = 0;
113  double x,y,z;
114  unsigned int node;
115  for (unsigned int i=0; i<Mesh->n_local_nodes_without_ghosts; ++i){
116  node = Mesh->local_nodes[i];
117  x = Mesh->nodes_coord[3*node+0];
118  y = Mesh->nodes_coord[3*node+1];
119  z = Mesh->nodes_coord[3*node+2];
120  if(x==0||y==0||z==0||x==1.0||y==1.0||z==1.0){
121  n_bc_dof+=3;
122  }
123  }
124 
125  int indbc = 0;
126  dof_on_boundary = new int [n_bc_dof];
127  for (unsigned int inode=0; inode<Mesh->n_local_nodes_without_ghosts; ++inode){
128  node = Mesh->local_nodes[inode];
129  x = Mesh->nodes_coord[3*node+0];
130  y = Mesh->nodes_coord[3*node+1];
131  z = Mesh->nodes_coord[3*node+2];
132  if(x==0||y==0||z==0||x==1.0||y==1.0||z==1.0){
133  dof_on_boundary[indbc+0] = 3*inode+0;
134  dof_on_boundary[indbc+1] = 3*inode+1;
135  dof_on_boundary[indbc+2] = 3*inode+2;
136  indbc+=3;
137  }
138  }
139  }
140 
141  void apply_dirichlet_conditions(Epetra_FECrsMatrix & K, Epetra_FEVector & F, double & displacement){
142  Epetra_MultiVector v(*StandardMap,true);
143  v.PutScalar(0.0);
144 
145  int node;
146  double x,y,z;
147  for (unsigned int inode=0; inode<Mesh->n_local_nodes_without_ghosts; ++inode){
148  node = Mesh->local_nodes[inode];
149  x = Mesh->nodes_coord[3*node+0];
150  y = Mesh->nodes_coord[3*node+1];
151  z = Mesh->nodes_coord[3*node+2];
152  if (x==0||y==0||z==0||x==1.0||y==1.0||z==1.0){
153  Epetra_SerialDenseVector u(3);
154  u = exactSolution(x,y,z);
155  v[0][StandardMap->LID(3*node+0)] = displacement*u(0);
156  v[0][StandardMap->LID(3*node+1)] = displacement*u(1);
157  v[0][StandardMap->LID(3*node+2)] = displacement*u(2);
158  }
159  }
160 
161  Epetra_MultiVector rhs_dir(*StandardMap,true);
162  K.Apply(v,rhs_dir);
163  F.Update(-1.0,rhs_dir,1.0);
164 
165  for (unsigned int inode=0; inode<Mesh->n_local_nodes_without_ghosts; ++inode){
166  node = Mesh->local_nodes[inode];
167  x = Mesh->nodes_coord[3*node+0];
168  y = Mesh->nodes_coord[3*node+1];
169  z = Mesh->nodes_coord[3*node+2];
170  if (x==0||y==0||z==0||x==1.0||y==1.0||z==1.0){
171  F[0][StandardMap->LID(3*node+0)] = v[0][StandardMap->LID(3*node+0)];
172  F[0][StandardMap->LID(3*node+1)] = v[0][StandardMap->LID(3*node+1)];
173  F[0][StandardMap->LID(3*node+2)] = v[0][StandardMap->LID(3*node+2)];
174  }
175  }
176  ML_Epetra::Apply_OAZToMatrix(dof_on_boundary,n_bc_dof,K);
177  }
178 
179  void get_material_parameters(unsigned int & e_lid, unsigned int & gp){
180  }
181 
182  void get_constitutive_tensors(Epetra_SerialDenseMatrix & deformation_gradient, Epetra_SerialDenseVector & piola_stress, Epetra_SerialDenseMatrix & tangent_piola){
183 
184  double det = deformation_gradient(0,0)*deformation_gradient(1,1)*deformation_gradient(2,2)
185  - deformation_gradient(0,0)*deformation_gradient(1,2)*deformation_gradient(2,1)
186  - deformation_gradient(0,1)*deformation_gradient(1,0)*deformation_gradient(2,2)
187  + deformation_gradient(0,1)*deformation_gradient(1,2)*deformation_gradient(2,0)
188  + deformation_gradient(0,2)*deformation_gradient(1,0)*deformation_gradient(2,1)
189  - deformation_gradient(0,2)*deformation_gradient(1,1)*deformation_gradient(2,0);
190 
191  Epetra_SerialDenseMatrix C(3,3);
192  Epetra_SerialDenseMatrix CC(3,3);
193  Epetra_SerialDenseMatrix ddJ5(6,6);
194  Epetra_SerialDenseVector LML(6);
195  Epetra_SerialDenseVector eye(6);
196  Epetra_SerialDenseVector dJ5(6);
197  Epetra_SerialDenseVector M(6);
198  Epetra_SerialDenseVector L(6);
199  Epetra_SerialDenseVector c(6);
200 
201  C.Multiply('T','N',1.0,deformation_gradient,deformation_gradient,0.0);
202  CC.Multiply('N','N',1.0,C,C,0.0);
203 
204  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);
205 
206  L(0) = (1.0/(det*det))*(C(1,1)*C(2,2)-C(1,2)*C(2,1));
207  L(1) = (1.0/(det*det))*(C(0,0)*C(2,2)-C(0,2)*C(2,0));
208  L(2) = (1.0/(det*det))*(C(0,0)*C(1,1)-C(0,1)*C(1,0));
209  L(3) = (1.0/(det*det))*(C(0,2)*C(1,0)-C(0,0)*C(1,2));
210  L(4) = (1.0/(det*det))*(C(0,1)*C(1,2)-C(0,2)*C(1,1));
211  L(5) = (1.0/(det*det))*(C(0,2)*C(2,1)-C(0,1)*C(2,2));
212 
213  M(0) = mu4*sin_plyagl*sin_plyagl+mu5*cos_plyagl*cos_plyagl;
214  M(1) = mu4*cos_plyagl*cos_plyagl+mu5*sin_plyagl*sin_plyagl;
215  M(2) = mu5;
216  M(3) = 0.0;
217  M(4) = 0.0;
218  M(5) = (mu4-mu5)*cos_plyagl*sin_plyagl;
219 
220  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);
221  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);
222  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);
223  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));
224  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));
225  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));
226 
227  eye(0) = 1.0; eye(1) = 1.0; eye(2) = 1.0; eye(3) = 0.0; eye(4) = 0.0; eye(5) = 0.0;
228 
229  double I1 = C(0,0) + C(1,1) + C(2,2);
230  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);
231  double I2 = (1.0/2.0)*(I1*I1-II1);
232  double I3 = det*det;
233  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);
234  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);
235  double J5 = I5 - I1*I4 + I2*trm;
236  double pI3 = std::pow(I3,-beta3);
237  double pI4 = std::pow(I4, beta4);
238  double pJ5 = std::pow(J5, beta5);
239 
240  for (unsigned int i=0; i<6; ++i){
241  dJ5(i) = J5*L(i) - I3*LML(i);
242  piola_stress(i) = 2.0*mu1*eye(i) + 2.0*mu2*(I1*eye(i)-c(i)) + (2.0*mu3*det*det-mu)*L(i)
243  + (2.0/ptrmbeta4)*pI4*M(i) + (2.0/ptrmbeta5)*pJ5*dJ5(i) - 2.0*trm*pI3*L(i);
244  }
245 
246  tensor_product(4.0*mu2,eye,eye,tangent_piola,0.0);
247  sym_tensor_product(-4.0*mu2,eye,eye,tangent_piola,1.0);
248 
249  tensor_product(4.0*mu3*det*det,L,L,tangent_piola,1.0);
250 
251  sym_tensor_product(-4.0*mu3*det*det+2.0*mu,L,L,tangent_piola,1.0);
252 
253  tensor_product(J5,L,L,ddJ5,0.0);
254  sym_tensor_product(-J5,L,L,ddJ5,1.0);
255 
256  tensor_product(-I3,L,LML,ddJ5,1.0);
257  tensor_product(-I3,LML,L,ddJ5,1.0);
258 
259  sym_tensor_product(I3,L,LML,ddJ5,1.0);
260  sym_tensor_product(I3,LML,L,ddJ5,1.0);
261 
262  tensor_product(4.0*beta4*pI4/(I4*ptrmbeta4),M,M,tangent_piola,1.0);
263  tensor_product(4.0*beta5*pJ5/(J5*ptrmbeta5),dJ5,dJ5,tangent_piola,1.0);
264 
265  tensor_product(4.0*trm*beta3*pI3,L,L,tangent_piola,1.0);
266  sym_tensor_product(4.0*trm*pI3,L,L,tangent_piola,1.0);
267 
268  ddJ5.Scale(4.0*pJ5/ptrmbeta5);
269  tangent_piola += ddJ5;
270  }
271 
272  Epetra_SerialDenseVector get_neumannBc(Epetra_SerialDenseMatrix & matrix_X, Epetra_SerialDenseMatrix & xg, unsigned int & gp){
273  Epetra_SerialDenseVector h(3);
274  std::cout << "**Err: Not using that method in this example!\n";
275  return h;
276  }
277 
278  Epetra_SerialDenseVector get_forcing(double & x1, double & x2, double & x3, unsigned int & e_lid, unsigned int & gp){
279  Epetra_SerialDenseVector f(3);
280  std::cout << "**Err: Not using that method in this example!\n";
281  return f;
282  }
283 
284  void get_material_parameters_for_recover(unsigned int & e_lid){
285  }
286 
287  void get_stress_for_recover(Epetra_SerialDenseMatrix & deformation_gradient, double & det, Epetra_SerialDenseMatrix & piola_stress){
288  }
289 
290 };
291 #endif
Epetra_Comm * Comm
u
Definition: run.m:9
void get_constitutive_tensors(Epetra_SerialDenseMatrix &deformation_gradient, Epetra_SerialDenseVector &piola_stress, Epetra_SerialDenseMatrix &tangent_piola)
std::vector< int > local_nodes
Definition: meshpp.hpp:80
void tensor_product(double scalarAB, Epetra_SerialDenseVector &A, Epetra_SerialDenseVector &B, Epetra_SerialDenseMatrix &AoB, double scalarThis)
void get_stress_for_recover(Epetra_SerialDenseMatrix &deformation_gradient, double &det, Epetra_SerialDenseMatrix &piola_stress)
std::vector< double > nodes_coord
Definition: meshpp.hpp:77
unsigned int n_bc_dof
Epetra_SerialDenseVector N_cells
Definition: meshpp.hpp:73
int n_local_cells
Definition: meshpp.hpp:94
Epetra_SerialDenseVector exactSolution(double &x1, double &x2, double &x3)
Epetra_SerialDenseVector detJac_cells
Definition: meshpp.hpp:73
void get_matrix_and_rhs(Epetra_Vector &x, Epetra_FECrsMatrix &K, Epetra_FEVector &F)
std::vector< int > local_cells
Definition: meshpp.hpp:81
clc clear all close all mesh
Definition: run.m:5
void assemblePureDirichlet_homogeneousForcing(Epetra_Vector &x, Epetra_FECrsMatrix &K, Epetra_FEVector &F)
int n_local_nodes_without_ghosts
Definition: meshpp.hpp:92
int el_type
Definition: meshpp.hpp:90
std::vector< int > cells_nodes
Definition: meshpp.hpp:78
void set_parameters(Epetra_SerialDenseVector &x, double &angle)
L
Definition: costFunction.m:66
Epetra_Import * ImportToOverlapMap
Epetra_Map * StandardMap
void apply_dirichlet_conditions(Epetra_FECrsMatrix &K, Epetra_FEVector &F, double &displacement)
unsigned int n_gauss_cells
Definition: meshpp.hpp:98
Epetra_SerialDenseVector get_forcing(double &x1, double &x2, double &x3, unsigned int &e_lid, unsigned int &gp)
void sym_tensor_product(double scalarAB, Epetra_SerialDenseVector &A, Epetra_SerialDenseVector &B, Epetra_SerialDenseMatrix &AoB, double scalarThis)
Epetra_SerialDenseVector get_neumannBc(Epetra_SerialDenseMatrix &matrix_X, Epetra_SerialDenseMatrix &xg, unsigned int &gp)
for i
Definition: costFunction.m:38
Epetra_Map * OverlapMap
compressibleHyperelasticity_linearPatchTest(Epetra_Comm &comm, Teuchos::ParameterList &Parameters)
Epetra_SerialDenseVector gauss_weight_cells
Definition: meshpp.hpp:99