Trilinos based (stochastic) FEM solvers
neumannInnerSurface_PolyconvexHGO.hpp
Go to the documentation of this file.
1 /*
2 Brian Staber (brian.staber@gmail.com)
3 */
4 
5 #ifndef NEUMANNINNERSURFACE_POLYCONVEXHGO_HPP
6 #define NEUMANNINNERSURFACE_POLYCONVEXHGO_HPP
7 
8 #include "tensor_calculus.hpp"
9 #include "laplacepp.hpp"
10 #include "nearlyIncompressibleHyperelasticity.hpp"
11 
12 class neumannInnerSurface_PolyconvexHGO : public nearlyIncompressibleHyperelasticity
13 {
14 public:
15 
16  laplace * Laplace;
17 
18  double mu1, mu2, mu3, mu4, beta3, beta4, theta;
19  Epetra_SerialDenseVector a, b;
20 
21  neumannInnerSurface_PolyconvexHGO(Epetra_Comm & comm, Teuchos::ParameterList & Parameters){
22 
23  std::string mesh_file = Teuchos::getParameter<std::string>(Parameters.sublist("Mesh"), "mesh_file");
24  std::string boundary_file = Teuchos::getParameter<std::string>(Parameters.sublist("Mesh"), "boundary_file");
25  unsigned int number_physical_groups = Teuchos::getParameter<unsigned int>(Parameters.sublist("Mesh"), "nb_phys_groups");
26  std::string select_model = Teuchos::getParameter<std::string>(Parameters.sublist("Mesh"), "model");
27 
28  mu1 = Teuchos::getParameter<double>(Parameters.sublist(select_model), "mu1");
29  mu2 = Teuchos::getParameter<double>(Parameters.sublist(select_model), "mu2");
30  mu3 = Teuchos::getParameter<double>(Parameters.sublist(select_model), "mu3");
31  mu4 = Teuchos::getParameter<double>(Parameters.sublist(select_model), "mu4");
32  beta3 = Teuchos::getParameter<double>(Parameters.sublist(select_model), "beta3");
33  beta4 = Teuchos::getParameter<double>(Parameters.sublist(select_model), "beta4");
34  theta = Teuchos::getParameter<double>(Parameters.sublist(select_model), "theta");
35 
36  Mesh = new mesh(comm, mesh_file, 1000.0);
37  Mesh->read_boundary_file(boundary_file,number_physical_groups);
38  Comm = Mesh->Comm;
39 
40  Laplace = new laplace(*Mesh);
41 
42  //$$Laplace
43  Epetra_FECrsMatrix matrix(Copy,*Laplace->FEGraph);
44  Epetra_Vector psi(*Laplace->StandardMap);
45  Epetra_Vector phi(*Laplace->StandardMap);
46  Epetra_FEVector rhs(*Laplace->StandardMap);
47  int bc_indx[2];
48  double bc_val[2];
49  bc_val[0] = 0.0; bc_val[1] = 1.0;
50  //Problem number one with aztec
51  bc_indx[0] = 2; bc_indx[1] = 3;
52  Laplace->solve_aztec(Parameters.sublist("Laplace"), matrix, phi, rhs, &bc_indx[0], &bc_val[0]);
53  Laplace->print_solution(phi, "laplace_inlet_to_outlet_aztec.mtx");
54  //Problem number one with aztec
55  bc_indx[0] = 0; bc_indx[1] = 1;
56  Laplace->solve_aztec(Parameters.sublist("Laplace"), matrix, psi, rhs, &bc_indx[0], &bc_val[0]);
57  Laplace->print_solution(psi, "laplace_inner_to_outer_aztec.mtx");
58  //Get local directions
59  Laplace->compute_local_directions(phi, psi);
60  Laplace->compute_center_local_directions(phi, psi);
61  //$$End
62 
63  StandardMap = new Epetra_Map(-1,3*Mesh->n_local_nodes_without_ghosts,&Mesh->local_dof_without_ghosts[0],0,*Comm);
64  OverlapMap = new Epetra_Map(-1,3*Mesh->n_local_nodes,&Mesh->local_dof[0],0,*Comm);
65  ImportToOverlapMap = new Epetra_Import(*OverlapMap,*StandardMap);
66  create_FECrsGraph();
67 
68  a.Resize(3);
69  b.Resize(3);
70  setup_dirichlet_conditions();
71  }
72 
73  ~neumannInnerSurface_PolyconvexHGO(){
74  }
75 
76  void get_matrix_and_rhs(Epetra_Vector & x, Epetra_FECrsMatrix & K, Epetra_FEVector & F){
77  assembleMixedDirichletDeformationDependentNeumann_homogeneousForcing(x,K,F);
78  }
79 
80  void setup_dirichlet_conditions(){
81  //setup_clamp();
82  setup_semislipbc();
83  //setup_slipbc_and_clamp();
84  //setup_slipbc();
85  }
86 
87  void apply_dirichlet_conditions(Epetra_FECrsMatrix & K, Epetra_FEVector & F, double & displacement){
88  //apply_clamp(K,F,displacement);
89  apply_semiselipbc(K,F,displacement);
90  //apply_slipbc_and_clamp(K,F,displacement);
91  //apply_slipbc(K,F,displacement);
92  }
93 
94  void get_material_parameters(unsigned int & e_lid, unsigned int & gp){
95 
96  int n_gauss_points = Mesh->n_gauss_cells;
97  int e_gid = Mesh->local_cells[e_lid];
98 
99  for (int i=0; i<3; ++i){
100  a(i) = cos(theta)*Laplace->laplace_direction_one(n_gauss_points*e_lid+gp,i)
101  + sin(theta)*Laplace->laplace_direction_two_cross_one(n_gauss_points*e_lid+gp,i);
102  b(i) = cos(theta)*Laplace->laplace_direction_one(n_gauss_points*e_lid+gp,i)
103  - sin(theta)*Laplace->laplace_direction_two_cross_one(n_gauss_points*e_lid+gp,i);
104  }
105 
106  }
107 
108  void get_constitutive_tensors_static_condensation(Epetra_SerialDenseMatrix & deformation_gradient,
109  double & det, Epetra_SerialDenseVector & inverse_cauchy,
110  Epetra_SerialDenseVector & piola_isc,
111  Epetra_SerialDenseVector & piola_vol,
112  Epetra_SerialDenseMatrix & tangent_piola_isc,
113  Epetra_SerialDenseMatrix & tangent_piola_vol){
114 
115  model_C(deformation_gradient, det, inverse_cauchy, piola_isc, piola_vol, tangent_piola_isc, tangent_piola_vol);
116 
117  }
118 
119  void get_internal_pressure(double & theta, double & pressure, double & dpressure){
120  double ptheta = std::pow(theta,beta3);
121  pressure = mu3*beta3*( (ptheta/theta) - (1.0/(ptheta*theta)) );
122  dpressure = mu3*beta3*( (beta3-1.0)*(ptheta/(theta*theta)) + (beta3+1.0)/(ptheta*theta*theta) );
123  }
124 
125  void get_material_parameters_for_recover(unsigned int & e_lid){
126 
127  for (int i=0; i<3; ++i){
128  a(i) = cos(theta)*Laplace->laplace_direction_one_center(e_lid,i)
129  + sin(theta)*Laplace->laplace_direction_two_cross_one_center(e_lid,i);
130  b(i) = cos(theta)*Laplace->laplace_direction_one_center(e_lid,i)
131  - sin(theta)*Laplace->laplace_direction_two_cross_one_center(e_lid,i);
132  }
133 
134  }
135 
136  void get_stress_for_recover(Epetra_SerialDenseMatrix & deformation_gradient, double & det, Epetra_SerialDenseMatrix & piola_stress){
137 
138  det = deformation_gradient(0,0)*deformation_gradient(1,1)*deformation_gradient(2,2)
139  - deformation_gradient(0,0)*deformation_gradient(1,2)*deformation_gradient(2,1)
140  - deformation_gradient(0,1)*deformation_gradient(1,0)*deformation_gradient(2,2)
141  + deformation_gradient(0,1)*deformation_gradient(1,2)*deformation_gradient(2,0)
142  + deformation_gradient(0,2)*deformation_gradient(1,0)*deformation_gradient(2,1)
143  - deformation_gradient(0,2)*deformation_gradient(1,1)*deformation_gradient(2,0);
144 
145  double alpha = std::pow(det,-2.0/3.0);
146  double beta = 1.0/(det*det);
147  Epetra_SerialDenseMatrix eye(3,3);
148  Epetra_SerialDenseMatrix M1(3,3), M2(3,3);
149  Epetra_SerialDenseMatrix C(3,3), L(3,3);
150  Epetra_SerialDenseMatrix piola_ani1(3,3), piola_ani2(3,3);
151 
152  eye(0,0) = 1.0; eye(0,1) = 0.0; eye(0,2) = 0.0;
153  eye(1,0) = 0.0; eye(1,1) = 1.0; eye(1,2) = 0.0;
154  eye(2,0) = 0.0; eye(2,1) = 0.0; eye(2,2) = 1.0;
155 
156  M1.Multiply('N','T',1.0,a,a,0.0);
157  M2.Multiply('N','T',1.0,b,b,0.0);
158 
159  C.Multiply('T','N',1.0,deformation_gradient,deformation_gradient,0.0);
160 
161  L(0,0) = (1.0/(det*det))*(C(1,1)*C(2,2)-C(1,2)*C(2,1));
162  L(1,1) = (1.0/(det*det))*(C(0,0)*C(2,2)-C(0,2)*C(2,0));
163  L(2,2) = (1.0/(det*det))*(C(0,0)*C(1,1)-C(0,1)*C(1,0));
164  L(1,2) = (1.0/(det*det))*(C(0,2)*C(1,0)-C(0,0)*C(1,2));
165  L(0,2) = (1.0/(det*det))*(C(0,1)*C(1,2)-C(0,2)*C(1,1));
166  L(0,1) = (1.0/(det*det))*(C(0,2)*C(2,1)-C(0,1)*C(2,2));
167  L(2,1) = L(1,2); L(2,0) = L(0,2); L(1,0) = L(0,1);
168 
169  double I1 = C(0,0) + C(1,1) + C(2,2);
170  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);
171  double I2 = (1.0/2.0)*(I1*I1-II1);
172  double I4_1 = C(0,0)*M1(0,0) + C(1,1)*M1(1,1) + C(2,2)*M1(2,2) + 2.0*C(0,1)*M1(0,1) + 2.0*C(0,2)*M1(0,2) + 2.0*C(1,2)*M1(1,2);
173  double I4_2 = C(0,0)*M2(0,0) + C(1,1)*M2(1,1) + C(2,2)*M2(2,2) + 2.0*C(0,1)*M2(0,1) + 2.0*C(0,2)*M2(0,2) + 2.0*C(1,2)*M2(1,2);
174  double pI2 = std::sqrt(I2);
175 
176  double S4_1 = (I4_1-1.0)*(I4_1-1.0);
177  double S4_2 = (I4_2-1.0)*(I4_2-1.0);
178 
179  double ptheta = std::pow(det,beta3);
180  double pressure = mu3*beta3*( (ptheta/det) - (1.0/(ptheta*det)) );
181 
182  for (unsigned int i=0; i<3; ++i){
183  for (unsigned int j=0; j<3; ++j){
184  piola_stress(i,j) = 2.0*mu1*alpha*(eye(i,j)-(1.0/3.0)*L(i,j))
185  + mu2*beta*( 3.0*pI2*(I1*eye(i,j)-C(i,j)) - 2.0*I2*pI2*L(i,j) )
186  + det*pressure*L(i,j);
187  piola_ani1(i,j) = 4.0*mu4*(I4_1-1.0)*exp(beta4*S4_1)*M1(i,j);
188  piola_ani2(i,j) = 4.0*mu4*(I4_2-1.0)*exp(beta4*S4_2)*M2(i,j);
189  }
190  }
191 
192  if (I4_1>1.0){
193  piola_stress += piola_ani1;
194  }
195  if (I4_2>1.0){
196  piola_stress += piola_ani2;
197  }
198 
199  }
200 
201  void setup_clamp(){
202  n_bc_dof = 0;
203  for (unsigned int i=0; i<Mesh->n_local_nodes_without_ghosts; ++i){
204  if (Mesh->nodes_to_boundaries(i,2)==1 || Mesh->nodes_to_boundaries(i,3)==1){
205  n_bc_dof+=3;
206  }
207  }
208 
209  int indbc = 0;
210  dof_on_boundary = new int [n_bc_dof];
211  for (unsigned int inode=0; inode<Mesh->n_local_nodes_without_ghosts; ++inode){
212  if (Mesh->nodes_to_boundaries(inode,2)==1 || Mesh->nodes_to_boundaries(inode,3)==1){
213  dof_on_boundary[indbc+0] = 3*inode+0;
214  dof_on_boundary[indbc+1] = 3*inode+1;
215  dof_on_boundary[indbc+2] = 3*inode+2;
216  indbc+=3;
217  }
218  }
219  }
220 
221  void setup_slipbc(){
222  const int nodesblk_gid0 = 480-1;
223  const int nodesblk_gid1 = 538-1;
224  const int nodesblk_gid2 = 577-1;
225 
226  int node;
227  n_bc_dof = 0;
228  for (unsigned int i=0; i<Mesh->n_local_nodes_without_ghosts; ++i){
229  if (Mesh->nodes_to_boundaries(i,2)==1 || Mesh->nodes_to_boundaries(i,3)==1){
230  node = Mesh->local_nodes[i];
231  switch (node){
232  case nodesblk_gid0:
233  n_bc_dof+=2;
234  break;
235  case nodesblk_gid1:
236  n_bc_dof+=2;
237  break;
238  case nodesblk_gid2:
239  n_bc_dof+=2;
240  break;
241  default:
242  n_bc_dof+=1;
243  break;
244  };
245  }
246  }
247 
248  int indbc = 0;
249  dof_on_boundary = new int [n_bc_dof];
250  for (unsigned int inode=0; inode<Mesh->n_local_nodes_without_ghosts; ++inode){
251  if (Mesh->nodes_to_boundaries(inode,2)==1 || Mesh->nodes_to_boundaries(inode,3)==1){
252  node = Mesh->local_nodes[inode];
253  switch (node){
254  case nodesblk_gid0:
255  dof_on_boundary[indbc+0] = 3*inode+0;
256  dof_on_boundary[indbc+1] = 3*inode+1;
257  indbc+=2;
258  break;
259  case nodesblk_gid1:
260  dof_on_boundary[indbc+0] = 3*inode+0;
261  dof_on_boundary[indbc+1] = 3*inode+2;
262  indbc+=2;
263  break;
264  case nodesblk_gid2:
265  dof_on_boundary[indbc+0] = 3*inode+0;
266  dof_on_boundary[indbc+1] = 3*inode+2;
267  indbc+=2;
268  break;
269  default:
270  dof_on_boundary[indbc+0] = 3*inode+0;
271  indbc+=1;
272  break;
273  };
274  }
275  }
276  }
277 
278  void setup_semislipbc(){
279  n_bc_dof = 0;
280  for (unsigned int i=0; i<Mesh->n_local_nodes_without_ghosts; ++i){
281  if (Mesh->nodes_to_boundaries(i,2)==1){
282  n_bc_dof+=1;
283  }
284  if (Mesh->nodes_to_boundaries(i,3)==1){
285  n_bc_dof+=3;
286  }
287  }
288 
289  int indbc = 0;
290  dof_on_boundary = new int [n_bc_dof];
291  for (unsigned int inode=0; inode<Mesh->n_local_nodes_without_ghosts; ++inode){
292  if (Mesh->nodes_to_boundaries(inode,2)==1){
293  dof_on_boundary[indbc] = 3*inode+0;
294  indbc+=1;
295  }
296  if (Mesh->nodes_to_boundaries(inode,3)==1){
297  dof_on_boundary[indbc+0] = 3*inode+0;
298  dof_on_boundary[indbc+1] = 3*inode+1;
299  dof_on_boundary[indbc+2] = 3*inode+2;
300  indbc+=3;
301  }
302  }
303  }
304 
305  void setup_slipbc_and_clamp(){
306  Epetra_IntSerialDenseVector nodesblk_gid(3);
307  nodesblk_gid(0) = 480-1;
308  nodesblk_gid(1) = 481-1;
309  nodesblk_gid(2) = 479-1;
310 
311  int node;
312  n_bc_dof = 0;
313  for (unsigned int i=0; i<Mesh->n_local_nodes_without_ghosts; ++i){
314  if (Mesh->nodes_to_boundaries(i,2)==1 || Mesh->nodes_to_boundaries(i,3)==1){
315  node = Mesh->local_nodes[i];
316  if (node==nodesblk_gid(0)||node==nodesblk_gid(1)||node==nodesblk_gid(2)){
317  n_bc_dof+=3;
318  }
319  else{
320  n_bc_dof+=1;
321  }
322  }
323  }
324 
325  int indbc = 0;
326  dof_on_boundary = new int [n_bc_dof];
327  for (unsigned int inode=0; inode<Mesh->n_local_nodes_without_ghosts; ++inode){
328  if (Mesh->nodes_to_boundaries(inode,2)==1 || Mesh->nodes_to_boundaries(inode,3)==1){
329  node = Mesh->local_nodes[inode];
330  if (node==nodesblk_gid(0)||node==nodesblk_gid(1)||node==nodesblk_gid(2)){
331  dof_on_boundary[indbc+0] = 3*inode+0;
332  dof_on_boundary[indbc+1] = 3*inode+1;
333  dof_on_boundary[indbc+2] = 3*inode+2;
334  indbc+=3;
335  }
336  else{
337  dof_on_boundary[indbc+0] = 3*inode+0;
338  indbc+=1;
339  }
340  }
341  }
342  }
343 
344  void apply_clamp(Epetra_FECrsMatrix & K, Epetra_FEVector & F, double & displacement){
345  if (n_bc_dof>0){
346  int node;
347  for (int inode=0; inode<Mesh->n_local_nodes_without_ghosts; ++inode){
348  node = Mesh->local_nodes[inode];
349  if (Mesh->nodes_to_boundaries(inode,2)==1 || Mesh->nodes_to_boundaries(inode,3)==1){
350  F[0][StandardMap->LID(3*node+0)] = displacement;
351  F[0][StandardMap->LID(3*node+1)] = displacement;
352  F[0][StandardMap->LID(3*node+2)] = displacement;
353  }
354  }
355  }
356  ML_Epetra::Apply_OAZToMatrix(dof_on_boundary,n_bc_dof,K);
357  }
358 
359  void apply_slipbc(Epetra_FECrsMatrix & K, Epetra_FEVector & F, double & displacement){
360  const int nodesblk_gid0 = 480-1;
361  const int nodesblk_gid1 = 538-1;
362  const int nodesblk_gid2 = 577-1;
363 
364  if (n_bc_dof>0){
365  int node;
366  for (int inode=0; inode<Mesh->n_local_nodes_without_ghosts; ++inode){
367  node = Mesh->local_nodes[inode];
368  if (Mesh->nodes_to_boundaries(inode,2)==1 || Mesh->nodes_to_boundaries(inode,3)==1){
369  switch (node){
370  case nodesblk_gid0:
371  F[0][StandardMap->LID(3*node+0)] = displacement;
372  F[0][StandardMap->LID(3*node+1)] = displacement;
373  break;
374  case nodesblk_gid1:
375  F[0][StandardMap->LID(3*node+0)] = displacement;
376  F[0][StandardMap->LID(3*node+2)] = displacement;
377  break;
378  case nodesblk_gid2:
379  F[0][StandardMap->LID(3*node+0)] = displacement;
380  F[0][StandardMap->LID(3*node+2)] = displacement;
381  break;
382  default:
383  F[0][StandardMap->LID(3*node+0)] = displacement;
384  break;
385  };
386  }
387  }
388  }
389  ML_Epetra::Apply_OAZToMatrix(dof_on_boundary,n_bc_dof,K);
390  }
391 
392  void apply_semiselipbc(Epetra_FECrsMatrix & K, Epetra_FEVector & F, double & displacement){
393  if (n_bc_dof>0){
394  int node;
395  for (int inode=0; inode<Mesh->n_local_nodes_without_ghosts; ++inode){
396  node = Mesh->local_nodes[inode];
397  if (Mesh->nodes_to_boundaries(inode,2)==1){
398  F[0][StandardMap->LID(3*node+0)] = displacement;
399  }
400  if (Mesh->nodes_to_boundaries(inode,3)==1){
401  F[0][StandardMap->LID(3*node+0)] = displacement;
402  F[0][StandardMap->LID(3*node+1)] = displacement;
403  F[0][StandardMap->LID(3*node+2)] = displacement;
404  }
405  }
406  }
407  ML_Epetra::Apply_OAZToMatrix(dof_on_boundary,n_bc_dof,K);
408  }
409 
410  void apply_slipbc_and_clamp(Epetra_FECrsMatrix & K, Epetra_FEVector & F, double & displacement){
411  Epetra_IntSerialDenseVector nodesblk_gid(3);
412  nodesblk_gid(0) = 480-1;
413  nodesblk_gid(1) = 481-1;
414  nodesblk_gid(2) = 479-1;
415 
416  if (n_bc_dof>0){
417  int node;
418  for (int inode=0; inode<Mesh->n_local_nodes_without_ghosts; ++inode){
419  node = Mesh->local_nodes[inode];
420  if (Mesh->nodes_to_boundaries(inode,2)==1 || Mesh->nodes_to_boundaries(inode,3)==1){
421  if (node==nodesblk_gid(0)||node==nodesblk_gid(1)||node==nodesblk_gid(2)){
422  F[0][StandardMap->LID(3*node+0)] = displacement;
423  F[0][StandardMap->LID(3*node+1)] = displacement;
424  F[0][StandardMap->LID(3*node+2)] = displacement;
425  }
426  else{
427  F[0][StandardMap->LID(3*node+0)] = displacement;
428  }
429  }
430  }
431  }
432  ML_Epetra::Apply_OAZToMatrix(dof_on_boundary,n_bc_dof,K);
433  }
434 
435  void model_C(Epetra_SerialDenseMatrix & deformation_gradient, double & det, Epetra_SerialDenseVector & L, Epetra_SerialDenseVector & piola_isc, Epetra_SerialDenseVector & piola_vol, Epetra_SerialDenseMatrix & tangent_piola_isc, Epetra_SerialDenseMatrix & tangent_piola_vol){
436 
437  det = deformation_gradient(0,0)*deformation_gradient(1,1)*deformation_gradient(2,2)
438  - deformation_gradient(0,0)*deformation_gradient(1,2)*deformation_gradient(2,1)
439  - deformation_gradient(0,1)*deformation_gradient(1,0)*deformation_gradient(2,2)
440  + deformation_gradient(0,1)*deformation_gradient(1,2)*deformation_gradient(2,0)
441  + deformation_gradient(0,2)*deformation_gradient(1,0)*deformation_gradient(2,1)
442  - deformation_gradient(0,2)*deformation_gradient(1,1)*deformation_gradient(2,0);
443 
444  double alpha = std::pow(det,-2.0/3.0);
445 
446  Epetra_SerialDenseVector eye(6);
447  Epetra_SerialDenseMatrix rightCauchy(3,3);
448  Epetra_SerialDenseVector M1(6), M2(6);
449  Epetra_SerialDenseVector C(6), D(6);
450  Epetra_SerialDenseVector piola_nh(6), piola_ani1(6), piola_ani2(6);
451 
452  M1(0) = a(0)*a(0); M2(0) = b(0)*b(0);
453  M1(1) = a(1)*a(1); M2(1) = b(1)*b(1);
454  M1(2) = a(2)*a(2); M2(2) = b(2)*b(2);
455  M1(3) = a(1)*a(2); M2(3) = b(1)*b(2);
456  M1(4) = a(0)*a(2); M2(4) = b(0)*b(2);
457  M1(5) = a(0)*a(1); M2(5) = b(0)*b(1);
458 
459  rightCauchy.Multiply('T','N',1.0,deformation_gradient,deformation_gradient,0.0);
460 
461  eye(0) = 1.0; eye(1) = 1.0; eye(2) = 1.0; eye(3) = 0.0; eye(4) = 0.0; eye(5) = 0.0;
462 
463  C(0) = rightCauchy(0,0); C(1) = rightCauchy(1,1); C(2) = rightCauchy(2,2);
464  C(3) = rightCauchy(1,2); C(4) = rightCauchy(0,2); C(5) = rightCauchy(0,1);
465 
466  L(0) = (1.0/(det*det))*(rightCauchy(1,1)*rightCauchy(2,2)-rightCauchy(1,2)*rightCauchy(2,1));
467  L(1) = (1.0/(det*det))*(rightCauchy(0,0)*rightCauchy(2,2)-rightCauchy(0,2)*rightCauchy(2,0));
468  L(2) = (1.0/(det*det))*(rightCauchy(0,0)*rightCauchy(1,1)-rightCauchy(0,1)*rightCauchy(1,0));
469  L(3) = (1.0/(det*det))*(rightCauchy(0,2)*rightCauchy(1,0)-rightCauchy(0,0)*rightCauchy(1,2));
470  L(4) = (1.0/(det*det))*(rightCauchy(0,1)*rightCauchy(1,2)-rightCauchy(0,2)*rightCauchy(1,1));
471  L(5) = (1.0/(det*det))*(rightCauchy(0,2)*rightCauchy(2,1)-rightCauchy(0,1)*rightCauchy(2,2));
472 
473  double I1 = C(0) + C(1) + C(2);
474  double II1 = C(0)*C(0) + C(1)*C(1) + C(2)*C(2) + 2.0*C(3)*C(3) + 2.0*C(4)*C(4) + 2.0*C(5)*C(5);
475  double I2 = (1.0/2.0)*(I1*I1-II1);
476  double I4_1 = C(0)*M1(0) + C(1)*M1(1) + C(2)*M1(2) + 2.0*C(5)*M1(5) + 2.0*C(4)*M1(4) + 2.0*C(3)*M1(3);
477  double I4_2 = C(0)*M2(0) + C(1)*M2(1) + C(2)*M2(2) + 2.0*C(5)*M2(5) + 2.0*C(4)*M2(4) + 2.0*C(3)*M2(3);
478  double pI2 = std::sqrt(I2);
479 
480  double S4_1 = (I4_1-1.0)*(I4_1-1.0);
481  double S4_2 = (I4_2-1.0)*(I4_2-1.0);
482 
483  for (unsigned int i=0; i<6; ++i){
484  D(i) = I1*eye(i) - C(i);
485  piola_nh(i) = 2.0*alpha*mu1*(eye(i)-(1.0/3.0)*I1*L(i));
486  piola_isc(i) = piola_nh(i) + (mu2/(det*det))*( 3.0*pI2*D(i) - 2.0*pI2*I2*L(i) );
487 
488  piola_ani1(i) = 4.0*mu4*(I4_1-1.0)*exp(beta4*S4_1)*M1(i);
489  piola_ani2(i) = 4.0*mu4*(I4_2-1.0)*exp(beta4*S4_2)*M2(i);
490 
491  piola_vol(i) = det*L(i);
492  }
493 
494  double scalarAB;
495 
496  //scalarAB = det;
497  tensor_product(det,L,L,tangent_piola_vol,0.0);
498  //scalarAB = -2.0*det;
499  sym_tensor_product(-2.0*det,L,L,tangent_piola_vol,1.0);
500 
501  //scalarAB = -2.0/3.0;
502  tensor_product(-2.0/3.0,piola_nh,L,tangent_piola_isc,0.0);
503  tensor_product(-2.0/3.0,L,piola_nh,tangent_piola_isc,1.0);
504 
505  //scalarAB = -6.0*mu2*pI2/(det*det);
506  tensor_product(-6.0*mu2*pI2/(det*det),D,eye,tangent_piola_isc,1.0);
507  tensor_product(-6.0*mu2*pI2/(det*det),eye,D,tangent_piola_isc,1.0);
508 
509  //scalarAB = (-4.0/9.0)*mu1*alpha*I1 + 4.0*mu2*pI2*I2/(det*det);
510  tensor_product((-4.0/9.0)*mu1*alpha*I1 + 4.0*mu2*pI2*I2/(det*det),L,L,tangent_piola_isc,1.0);
511 
512  //scalarAB = (4.0/3.0)*mu1*alpha*I1 + 4.0*mu2*pI2*I2/(det*det);
513  sym_tensor_product((4.0/3.0)*mu1*alpha*I1 + 4.0*mu2*pI2*I2/(det*det),L,L,tangent_piola_isc,1.0);
514 
515  //scalarAB = 3.0*mu2/(det*det*pI2);
516  tensor_product(3.0*mu2/(det*det*pI2),D,D,tangent_piola_isc,1.0);
517 
518  //scalarAB = 6.0*mu2*pI2/(det*det);
519  tensor_product(6.0*mu2*pI2/(det*det),eye,eye,tangent_piola_isc,1.0);
520 
521  //scalarAB = -scalarAB;
522  sym_tensor_product(-6.0*mu2*pI2/(det*det),eye,eye,tangent_piola_isc,1.0);
523 
524  if (I4_1>1.0){
525  piola_isc += piola_ani1;
526  //scalarAB = (8.0*mu4 + 16.0*mu4*beta4*S4_1)*exp(beta4*S4_1);
527  tensor_product((8.0*mu4 + 16.0*mu4*beta4*S4_1)*exp(beta4*S4_1),M1,M1,tangent_piola_isc,1.0);
528  }
529  if (I4_2>1.0){
530  piola_isc += piola_ani2;
531  //scalarAB = (8.0*mu4 + 16.0*mu4*beta4*S4_2)*exp(beta4*S4_2);
532  tensor_product((8.0*mu4 + 16.0*mu4*beta4*S4_2)*exp(beta4*S4_2),M2,M2,tangent_piola_isc,1.0);
533  }
534  }
535 
536  void model_B(Epetra_SerialDenseMatrix & deformation_gradient, double & det, Epetra_SerialDenseVector & piola_isc, Epetra_SerialDenseVector & piola_vol, Epetra_SerialDenseMatrix & tangent_piola_isc, Epetra_SerialDenseMatrix & tangent_piola_vol){
537 
538  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);
539 
540  double alpha = std::pow(det,-2.0/3.0);
541 
542  Epetra_SerialDenseMatrix rightCauchy(3,3);
543  Epetra_SerialDenseVector M1(6), M2(6);
544  Epetra_SerialDenseVector eye(6);
545  Epetra_SerialDenseVector C(6), D(6), L(6);
546  Epetra_SerialDenseVector dI4_1(6), dI4_2(6);
547  Epetra_SerialDenseVector piola_nh(6), piola_ani1(6), piola_ani2(6);
548 
549  M1(0) = a(0)*a(0); M2(0) = b(0)*b(0);
550  M1(1) = a(1)*a(1); M2(1) = b(1)*b(1);
551  M1(2) = a(2)*a(2); M2(2) = b(2)*b(2);
552  M1(3) = a(1)*a(2); M2(3) = b(1)*b(2);
553  M1(4) = a(0)*a(2); M2(4) = b(0)*b(2);
554  M1(5) = a(0)*a(1); M2(5) = b(0)*b(1);
555 
556  rightCauchy.Multiply('T','N',1.0,deformation_gradient,deformation_gradient,0.0);
557 
558  eye(0) = 1.0; eye(1) = 1.0; eye(2) = 1.0; eye(3) = 0.0; eye(4) = 0.0; eye(5) = 0.0;
559 
560  C(0) = rightCauchy(0,0); C(1) = rightCauchy(1,1); C(2) = rightCauchy(2,2);
561  C(3) = rightCauchy(1,2); C(4) = rightCauchy(0,2); C(5) = rightCauchy(0,1);
562 
563  L(0) = (1.0/(det*det))*(rightCauchy(1,1)*rightCauchy(2,2)-rightCauchy(1,2)*rightCauchy(2,1));
564  L(1) = (1.0/(det*det))*(rightCauchy(0,0)*rightCauchy(2,2)-rightCauchy(0,2)*rightCauchy(2,0));
565  L(2) = (1.0/(det*det))*(rightCauchy(0,0)*rightCauchy(1,1)-rightCauchy(0,1)*rightCauchy(1,0));
566  L(3) = (1.0/(det*det))*(rightCauchy(0,2)*rightCauchy(1,0)-rightCauchy(0,0)*rightCauchy(1,2));
567  L(4) = (1.0/(det*det))*(rightCauchy(0,1)*rightCauchy(1,2)-rightCauchy(0,2)*rightCauchy(1,1));
568  L(5) = (1.0/(det*det))*(rightCauchy(0,2)*rightCauchy(2,1)-rightCauchy(0,1)*rightCauchy(2,2));
569 
570  double I1 = C(0) + C(1) + C(2);
571  double II1 = C(0)*C(0) + C(1)*C(1) + C(2)*C(2) + 2.0*C(3)*C(3) + 2.0*C(4)*C(4) + 2.0*C(5)*C(5);
572  double I2 = (1.0/2.0)*(I1*I1-II1);
573  double I4_1 = C(0)*M1(0) + C(1)*M1(1) + C(2)*M1(2) + 2.0*C(5)*M1(5) + 2.0*C(4)*M1(4) + 2.0*C(3)*M1(3);
574  double I4_2 = C(0)*M2(0) + C(1)*M2(1) + C(2)*M2(2) + 2.0*C(5)*M2(5) + 2.0*C(4)*M2(4) + 2.0*C(3)*M2(3);
575  double pI2 = std::sqrt(I2);
576 
577  double S4_1 = (alpha*I4_1-1.0)*(alpha*I4_1-1.0);
578  double S4_2 = (alpha*I4_2-1.0)*(alpha*I4_2-1.0);
579 
580  for (unsigned int i=0; i<6; ++i){
581  D(i) = I1*eye(i) - C(i);
582  dI4_1(i) = M1(i)-(1.0/3.0)*I4_1*L(i);
583  dI4_2(i) = M2(i)-(1.0/3.0)*I4_2*L(i);
584  piola_nh(i) = 2.0*alpha*mu1*(eye(i)-(1.0/3.0)*I1*L(i));
585  piola_isc(i) = piola_nh(i) + (mu2/(det*det))*( 3.0*pI2*D(i) - 2.0*pI2*I2*L(i) );
586 
587  piola_ani1(i) = 4.0*mu4*alpha*(alpha*I4_1-1.0)*exp(beta4*S4_1)*(M1(i)-(1.0/3.0)*I4_1*L(i));
588  piola_ani2(i) = 4.0*mu4*alpha*(alpha*I4_2-1.0)*exp(beta4*S4_2)*(M2(i)-(1.0/3.0)*I4_2*L(i));;
589 
590  piola_vol(i) = det*L(i);
591  }
592 
593  double scalarAB;
594 
595  scalarAB = det;
596  tensor_product(det,L,L,tangent_piola_vol,0.0);
597  scalarAB = -2.0*det;
598  sym_tensor_product(scalarAB,L,L,tangent_piola_vol,1.0);
599 
600  scalarAB = -2.0/3.0;
601  tensor_product(scalarAB,piola_nh,L,tangent_piola_isc,0.0);
602  tensor_product(scalarAB,L,piola_nh,tangent_piola_isc,1.0);
603 
604  scalarAB = -6.0*mu2*pI2/(det*det);
605  tensor_product(scalarAB,D,eye,tangent_piola_isc,1.0);
606  tensor_product(scalarAB,eye,D,tangent_piola_isc,1.0);
607 
608  scalarAB = (-4.0/9.0)*mu1*alpha*I1 + 4.0*mu2*pI2*I2/(det*det);
609  tensor_product(scalarAB,L,L,tangent_piola_isc,1.0);
610 
611  scalarAB = (4.0/3.0)*mu1*alpha*I1 + 4.0*mu2*pI2*I2/(det*det);
612  sym_tensor_product(scalarAB,L,L,tangent_piola_isc,1.0);
613 
614  scalarAB = 3.0*mu2/(det*det*pI2);
615  tensor_product(scalarAB,D,D,tangent_piola_isc,1.0);
616 
617  scalarAB = 6.0*mu2*pI2/(det*det);
618  tensor_product(scalarAB,eye,eye,tangent_piola_isc,1.0);
619 
620  scalarAB = -scalarAB;
621  sym_tensor_product(scalarAB,eye,eye,tangent_piola_isc,1.0);
622 
623  if (I4_1>1.0){
624  piola_isc += piola_ani1;
625  scalarAB = -(8.0/3.0)*mu4*alpha*(alpha*I4_1-1.0)*exp(beta4*S4_1);
626  tensor_product(scalarAB,dI4_1,L,tangent_piola_isc,1.0);
627  scalarAB = (8.0/3.0)*mu4*alpha*alpha*exp(beta4*S4_1);
628  tensor_product(scalarAB,dI4_1,dI4_1,tangent_piola_isc,1.0);
629  scalarAB = 16.0*mu4*beta4*alpha*alpha*S4_1*exp(beta4*S4_1);
630  tensor_product(scalarAB,dI4_1,dI4_1,tangent_piola_isc,1.0);
631  }
632  if (I4_2>1.0){
633  piola_isc += piola_ani2;
634  scalarAB = -(8.0/3.0)*mu4*alpha*(alpha*I4_2-1.0)*exp(beta4*S4_2);
635  tensor_product(scalarAB,dI4_2,L,tangent_piola_isc,1.0);
636  scalarAB = (8.0/3.0)*mu4*alpha*alpha*exp(beta4*S4_2);
637  tensor_product(scalarAB,dI4_2,dI4_2,tangent_piola_isc,1.0);
638  scalarAB = 16.0*mu4*beta4*alpha*alpha*S4_2*exp(beta4*S4_2);
639  tensor_product(scalarAB,dI4_2,dI4_2,tangent_piola_isc,1.0);
640  }
641 
642  }
643 
644  void model_A(Epetra_SerialDenseMatrix & deformation_gradient, double & det, Epetra_SerialDenseVector & piola_isc, Epetra_SerialDenseVector & piola_vol, Epetra_SerialDenseMatrix & tangent_piola_isc, Epetra_SerialDenseMatrix & tangent_piola_vol){
645 
646  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);
647 
648  double alpha = std::pow(det,-2.0/3.0);
649 
650  Epetra_SerialDenseMatrix rightCauchy(3,3);
651  Epetra_SerialDenseVector M1(6), M2(6);
652  Epetra_SerialDenseVector eye(6);
653  Epetra_SerialDenseVector C(6), L(6), CC(6);
654  Epetra_SerialDenseVector CMMC1(6), CMMC2(6);
655  Epetra_SerialDenseVector D(6);
656  Epetra_SerialDenseVector dK3_1(6), dK3_2(6);
657  Epetra_SerialDenseVector piola_nh(6), piola_ani1(6), piola_ani2(6);
658 
659  M1(0) = a(0)*a(0);
660  M1(1) = a(1)*a(1);
661  M1(2) = a(2)*a(2);
662  M1(3) = a(1)*a(2);
663  M1(4) = a(0)*a(2);
664  M1(5) = a(0)*a(1);
665 
666  M2(0) = b(0)*b(0);
667  M2(1) = b(1)*b(1);
668  M2(2) = b(2)*b(2);
669  M2(3) = b(1)*b(2);
670  M2(4) = b(0)*b(2);
671  M2(5) = b(0)*b(1);
672 
673  rightCauchy.Multiply('T','N',1.0,deformation_gradient,deformation_gradient,0.0);
674 
675  eye(0) = 1.0; eye(1) = 1.0; eye(2) = 1.0; eye(3) = 0.0; eye(4) = 0.0; eye(5) = 0.0;
676 
677  C(0) = rightCauchy(0,0); C(1) = rightCauchy(1,1); C(2) = rightCauchy(2,2);
678  C(3) = rightCauchy(1,2); C(4) = rightCauchy(0,2); C(5) = rightCauchy(0,1);
679 
680  L(0) = (1.0/(det*det))*(rightCauchy(1,1)*rightCauchy(2,2)-rightCauchy(1,2)*rightCauchy(2,1));
681  L(1) = (1.0/(det*det))*(rightCauchy(0,0)*rightCauchy(2,2)-rightCauchy(0,2)*rightCauchy(2,0));
682  L(2) = (1.0/(det*det))*(rightCauchy(0,0)*rightCauchy(1,1)-rightCauchy(0,1)*rightCauchy(1,0));
683  L(3) = (1.0/(det*det))*(rightCauchy(0,2)*rightCauchy(1,0)-rightCauchy(0,0)*rightCauchy(1,2));
684  L(4) = (1.0/(det*det))*(rightCauchy(0,1)*rightCauchy(1,2)-rightCauchy(0,2)*rightCauchy(1,1));
685  L(5) = (1.0/(det*det))*(rightCauchy(0,2)*rightCauchy(2,1)-rightCauchy(0,1)*rightCauchy(2,2));
686 
687  CMMC1(0) = 2.0*C(0)*M1(0) + 2.0*C(4)*M1(4) + 2.0*C(5)*M1(5);
688  CMMC1(1) = 2.0*C(1)*M1(1) + 2.0*C(3)*M1(3) + 2.0*C(5)*M1(5);
689  CMMC1(2) = 2.0*C(2)*M1(2) + 2.0*C(3)*M1(3) + 2.0*C(4)*M1(4);
690  CMMC1(3) = C(1)*M1(3) + C(3)*M1(1) + C(2)*M1(3) + C(3)*M1(2) + C(4)*M1(5) + C(5)*M1(4);
691  CMMC1(4) = C(0)*M1(4) + C(4)*M1(0) + C(2)*M1(4) + C(4)*M1(2) + C(3)*M1(5) + C(5)*M1(3);
692  CMMC1(5) = C(0)*M1(5) + C(5)*M1(0) + C(1)*M1(5) + C(5)*M1(1) + C(3)*M1(4) + C(4)*M1(3);
693 
694  CMMC2(0) = 2.0*C(0)*M2(0) + 2.0*C(4)*M2(4) + 2.0*C(5)*M2(5);
695  CMMC2(1) = 2.0*C(1)*M2(1) + 2.0*C(3)*M2(3) + 2.0*C(5)*M2(5);
696  CMMC2(2) = 2.0*C(2)*M2(2) + 2.0*C(3)*M2(3) + 2.0*C(4)*M2(4);
697  CMMC2(3) = C(1)*M2(3) + C(3)*M2(1) + C(2)*M2(3) + C(3)*M2(2) + C(4)*M2(5) + C(5)*M2(4);
698  CMMC2(4) = C(0)*M2(4) + C(4)*M2(0) + C(2)*M2(4) + C(4)*M2(2) + C(3)*M2(5) + C(5)*M2(3);
699  CMMC2(5) = C(0)*M2(5) + C(5)*M2(0) + C(1)*M2(5) + C(5)*M2(1) + C(3)*M2(4) + C(4)*M2(3);
700 
701  CC(0) = C(0)*C(0) + C(4)*C(4) + C(5)*C(5);
702  CC(1) = C(1)*C(1) + C(3)*C(3) + C(5)*C(5);
703  CC(2) = C(2)*C(2) + C(3)*C(3) + C(4)*C(4);
704  CC(3) = C(1)*C(3) + C(2)*C(3) + C(4)*C(5);
705  CC(4) = C(0)*C(4) + C(2)*C(4) + C(3)*C(5);
706  CC(5) = C(0)*C(5) + C(1)*C(5) + C(3)*C(4);
707 
708  double I1 = C(0) + C(1) + C(2);
709  double II1 = C(0)*C(0) + C(1)*C(1) + C(2)*C(2) + 2.0*C(3)*C(3) + 2.0*C(4)*C(4) + 2.0*C(5)*C(5);
710  double I2 = (1.0/2.0)*(I1*I1-II1);
711  double I4_1 = C(0)*M1(0) + C(1)*M1(1) + C(2)*M1(2) + 2.0*C(5)*M1(5) + 2.0*C(4)*M1(4) + 2.0*C(3)*M1(3);
712  double I4_2 = C(0)*M2(0) + C(1)*M2(1) + C(2)*M2(2) + 2.0*C(5)*M2(5) + 2.0*C(4)*M2(4) + 2.0*C(3)*M2(3);
713  double I5_1 = CC(0)*M1(0) + CC(1)*M1(1) + CC(2)*M1(2) + 2.0*CC(5)*M1(5) + 2.0*CC(4)*M1(4) + 2.0*CC(3)*M1(3);
714  double I5_2 = CC(0)*M2(0) + CC(1)*M2(1) + CC(2)*M2(2) + 2.0*CC(5)*M2(5) + 2.0*CC(4)*M2(4) + 2.0*CC(3)*M2(3);
715  double K3_1 = I1*I4_1-I5_1;
716  double K3_2 = I1*I4_2-I5_2;
717  double pI2 = std::sqrt(I2);
718 
719  for (unsigned int i=0; i<6; ++i){
720  D(i) = I1*eye(i) - C(i);
721  piola_nh(i) = 2.0*alpha*mu1*(eye(i)-(1.0/3.0)*I1*L(i));
722  piola_isc(i) = piola_nh(i) + (mu2/(det*det))*( 3.0*pI2*D(i) - 2.0*pI2*I2*L(i) );
723 
724  dK3_1(i) = I4_1*eye(i) + I1*M1(i) - CMMC1(i);
725  dK3_2(i) = I4_2*eye(i) + I1*M2(i) - CMMC2(i);
726  piola_ani1(i) = 2.0*mu4*beta4*std::pow(K3_1-2.0,beta4-1.0)*dK3_1(i);
727  piola_ani2(i) = 2.0*mu4*beta4*std::pow(K3_2-2.0,beta4-1.0)*dK3_2(i);
728 
729  piola_vol(i) = det*L(i);
730  }
731 
732  double scalarAB;
733 
734  scalarAB = det;
735  tensor_product(det,L,L,tangent_piola_vol,0.0);
736  scalarAB = -2.0*det;
737  sym_tensor_product(scalarAB,L,L,tangent_piola_vol,1.0);
738 
739  scalarAB = -2.0/3.0;
740  tensor_product(scalarAB,piola_nh,L,tangent_piola_isc,0.0);
741  tensor_product(scalarAB,L,piola_nh,tangent_piola_isc,1.0);
742 
743  scalarAB = -6.0*mu2*pI2/(det*det);
744  tensor_product(scalarAB,D,eye,tangent_piola_isc,1.0);
745  tensor_product(scalarAB,eye,D,tangent_piola_isc,1.0);
746 
747  scalarAB = (-4.0/9.0)*mu1*alpha*I1 + 4.0*mu2*pI2*I2/(det*det);
748  tensor_product(scalarAB,L,L,tangent_piola_isc,1.0);
749 
750  scalarAB = (4.0/3.0)*mu1*alpha*I1 + 4.0*mu2*pI2*I2/(det*det);
751  sym_tensor_product(scalarAB,L,L,tangent_piola_isc,1.0);
752 
753  scalarAB = 3.0*mu2/(det*det*pI2);
754  tensor_product(scalarAB,D,D,tangent_piola_isc,1.0);
755 
756  scalarAB = 6.0*mu2*pI2/(det*det);
757  tensor_product(scalarAB,eye,eye,tangent_piola_isc,1.0);
758 
759  scalarAB = -scalarAB;
760  sym_tensor_product(scalarAB,eye,eye,tangent_piola_isc,1.0);
761 
762  if (K3_1>2.0){
763  piola_isc += piola_ani1;
764  scalarAB = 4.0*mu4*beta4*(beta4-1.0)*std::pow(K3_1-2.0,beta4-2.0);
765  tensor_product(scalarAB,dK3_1,dK3_1,tangent_piola_isc,1.0);
766 
767  scalarAB = 4.0*mu4*beta4*std::pow(K3_1-2.0,beta4-1.0);
768  tensor_product(scalarAB,eye,M1,tangent_piola_isc,1.0);
769  tensor_product(scalarAB,M1,eye,tangent_piola_isc,1.0);
770  scalarAB = -scalarAB;
771  sym_tensor_product(scalarAB,M1,eye,tangent_piola_isc,1.0);
772  sym_tensor_product(scalarAB,eye,M1,tangent_piola_isc,1.0);
773  }
774  if (K3_2>2.0){
775  piola_isc += piola_ani2;
776  scalarAB = 4.0*mu4*beta4*(beta4-1.0)*std::pow(K3_2-2.0,beta4-2.0);
777  tensor_product(scalarAB,dK3_2,dK3_2,tangent_piola_isc,1.0);
778 
779  scalarAB = 4.0*mu4*beta4*std::pow(K3_2-2.0,beta4-1.0);
780  tensor_product(scalarAB,eye,M2,tangent_piola_isc,1.0);
781  tensor_product(scalarAB,M2,eye,tangent_piola_isc,1.0);
782  scalarAB = -scalarAB;
783  sym_tensor_product(scalarAB,M2,eye,tangent_piola_isc,1.0);
784  sym_tensor_product(scalarAB,eye,M2,tangent_piola_isc,1.0);
785  }
786  }
787 
788 };
789 
790 #endif
neumannInnerSurface_PolyconvexHGO::neumannInnerSurface_PolyconvexHGO
neumannInnerSurface_PolyconvexHGO(Epetra_Comm &comm, Teuchos::ParameterList &Parameters)
Definition: neumannInnerSurface_PolyconvexHGO.hpp:21
neumannInnerSurface_PolyconvexHGO::mu2
double mu2
Definition: neumannInnerSurface_PolyconvexHGO.hpp:18
neumannInnerSurface_PolyconvexHGO::model_C
void model_C(Epetra_SerialDenseMatrix &deformation_gradient, double &det, Epetra_SerialDenseVector &L, Epetra_SerialDenseVector &piola_isc, Epetra_SerialDenseVector &piola_vol, Epetra_SerialDenseMatrix &tangent_piola_isc, Epetra_SerialDenseMatrix &tangent_piola_vol)
Definition: neumannInnerSurface_PolyconvexHGO.hpp:435
tensor_calculus.hpp
baseClassFEM::Comm
Epetra_Comm * Comm
Definition: baseClassFEM.hpp:15
j
for j
Definition: costFunction.m:42
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
neumannInnerSurface_PolyconvexHGO::Laplace
laplace * Laplace
Definition: neumannInnerSurface_PolyconvexHGO.hpp:16
neumannInnerSurface_PolyconvexHGO::b
Epetra_SerialDenseVector b
Definition: neumannInnerSurface_PolyconvexHGO.hpp:19
neumannInnerSurface_PolyconvexHGO::apply_slipbc
void apply_slipbc(Epetra_FECrsMatrix &K, Epetra_FEVector &F, double &displacement)
Definition: neumannInnerSurface_PolyconvexHGO.hpp:359
neumannInnerSurface_PolyconvexHGO::setup_slipbc
void setup_slipbc()
Definition: neumannInnerSurface_PolyconvexHGO.hpp:221
hyperelasticity::n_bc_dof
unsigned int n_bc_dof
Definition: hyperelasticity.hpp:23
neumannInnerSurface_PolyconvexHGO::apply_semiselipbc
void apply_semiselipbc(Epetra_FECrsMatrix &K, Epetra_FEVector &F, double &displacement)
Definition: neumannInnerSurface_PolyconvexHGO.hpp:392
neumannInnerSurface_PolyconvexHGO::beta3
double beta3
Definition: neumannInnerSurface_PolyconvexHGO.hpp:18
laplace::print_solution
int print_solution(Epetra_Vector &lhs, std::string fileName)
Definition: laplacepp.cpp:399
laplace::compute_center_local_directions
void compute_center_local_directions(Epetra_Vector &laplace_one, Epetra_Vector &laplace_two)
Definition: laplacepp.cpp:307
neumannInnerSurface_PolyconvexHGO::get_material_parameters_for_recover
void get_material_parameters_for_recover(unsigned int &e_lid)
Definition: neumannInnerSurface_PolyconvexHGO.hpp:125
neumannInnerSurface_PolyconvexHGO::model_A
void model_A(Epetra_SerialDenseMatrix &deformation_gradient, double &det, Epetra_SerialDenseVector &piola_isc, Epetra_SerialDenseVector &piola_vol, Epetra_SerialDenseMatrix &tangent_piola_isc, Epetra_SerialDenseMatrix &tangent_piola_vol)
Definition: neumannInnerSurface_PolyconvexHGO.hpp:644
neumannInnerSurface_PolyconvexHGO::mu4
double mu4
Definition: neumannInnerSurface_PolyconvexHGO.hpp:18
neumannInnerSurface_PolyconvexHGO::a
Epetra_SerialDenseVector a
Definition: neumannInnerSurface_PolyconvexHGO.hpp:19
neumannInnerSurface_PolyconvexHGO::get_matrix_and_rhs
void get_matrix_and_rhs(Epetra_Vector &x, Epetra_FECrsMatrix &K, Epetra_FEVector &F)
Definition: neumannInnerSurface_PolyconvexHGO.hpp:76
neumannInnerSurface_PolyconvexHGO::mu3
double mu3
Definition: neumannInnerSurface_PolyconvexHGO.hpp:18
neumannInnerSurface_PolyconvexHGO::model_B
void model_B(Epetra_SerialDenseMatrix &deformation_gradient, double &det, Epetra_SerialDenseVector &piola_isc, Epetra_SerialDenseVector &piola_vol, Epetra_SerialDenseMatrix &tangent_piola_isc, Epetra_SerialDenseMatrix &tangent_piola_vol)
Definition: neumannInnerSurface_PolyconvexHGO.hpp:536
baseClassFEM::FEGraph
Epetra_FECrsGraph * FEGraph
Definition: baseClassFEM.hpp:20
neumannInnerSurface_PolyconvexHGO::get_internal_pressure
void get_internal_pressure(double &theta, double &pressure, double &dpressure)
Definition: neumannInnerSurface_PolyconvexHGO.hpp:119
mesh::nodes_to_boundaries
Epetra_IntSerialDenseMatrix nodes_to_boundaries
Definition: meshpp.hpp:75
mesh::local_cells
std::vector< int > local_cells
Definition: meshpp.hpp:81
mesh
clc clear all close all mesh
Definition: run.m:5
neumannInnerSurface_PolyconvexHGO::~neumannInnerSurface_PolyconvexHGO
~neumannInnerSurface_PolyconvexHGO()
Definition: neumannInnerSurface_PolyconvexHGO.hpp:73
neumannInnerSurface_PolyconvexHGO::setup_semislipbc
void setup_semislipbc()
Definition: neumannInnerSurface_PolyconvexHGO.hpp:278
mesh::n_local_nodes_without_ghosts
int n_local_nodes_without_ghosts
Definition: meshpp.hpp:92
neumannInnerSurface_PolyconvexHGO::get_constitutive_tensors_static_condensation
void get_constitutive_tensors_static_condensation(Epetra_SerialDenseMatrix &deformation_gradient, double &det, Epetra_SerialDenseVector &inverse_cauchy, Epetra_SerialDenseVector &piola_isc, Epetra_SerialDenseVector &piola_vol, Epetra_SerialDenseMatrix &tangent_piola_isc, Epetra_SerialDenseMatrix &tangent_piola_vol)
Definition: neumannInnerSurface_PolyconvexHGO.hpp:108
L
L
Definition: costFunction.m:66
laplace::laplace_direction_one
Epetra_SerialDenseMatrix laplace_direction_one
Definition: laplacepp.hpp:28
baseClassFEM::ImportToOverlapMap
Epetra_Import * ImportToOverlapMap
Definition: baseClassFEM.hpp:19
neumannInnerSurface_PolyconvexHGO::get_material_parameters
void get_material_parameters(unsigned int &e_lid, unsigned int &gp)
Definition: neumannInnerSurface_PolyconvexHGO.hpp:94
baseClassFEM::StandardMap
Epetra_Map * StandardMap
Definition: baseClassFEM.hpp:18
mesh::n_gauss_cells
unsigned int n_gauss_cells
Definition: meshpp.hpp:98
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
laplace::laplace_direction_two_cross_one
Epetra_SerialDenseMatrix laplace_direction_two_cross_one
Definition: laplacepp.hpp:30
neumannInnerSurface_PolyconvexHGO::theta
double theta
Definition: neumannInnerSurface_PolyconvexHGO.hpp:18
laplace::laplace_direction_one_center
Epetra_SerialDenseMatrix laplace_direction_one_center
Definition: laplacepp.hpp:28
neumannInnerSurface_PolyconvexHGO::apply_clamp
void apply_clamp(Epetra_FECrsMatrix &K, Epetra_FEVector &F, double &displacement)
Definition: neumannInnerSurface_PolyconvexHGO.hpp:344
neumannInnerSurface_PolyconvexHGO::get_stress_for_recover
void get_stress_for_recover(Epetra_SerialDenseMatrix &deformation_gradient, double &det, Epetra_SerialDenseMatrix &piola_stress)
Definition: neumannInnerSurface_PolyconvexHGO.hpp:136
neumannInnerSurface_PolyconvexHGO
Definition: neumannInnerSurface_PolyconvexHGO.hpp:12
laplace::laplace_direction_two_cross_one_center
Epetra_SerialDenseMatrix laplace_direction_two_cross_one_center
Definition: laplacepp.hpp:30
neumannInnerSurface_PolyconvexHGO::setup_slipbc_and_clamp
void setup_slipbc_and_clamp()
Definition: neumannInnerSurface_PolyconvexHGO.hpp:305
laplace::compute_local_directions
void compute_local_directions(Epetra_Vector &laplace_one, Epetra_Vector &laplace_two)
Definition: laplacepp.cpp:229
neumannInnerSurface_PolyconvexHGO::apply_slipbc_and_clamp
void apply_slipbc_and_clamp(Epetra_FECrsMatrix &K, Epetra_FEVector &F, double &displacement)
Definition: neumannInnerSurface_PolyconvexHGO.hpp:410
beta
modelParameters beta
Definition: run_station15.m:10
nearlyIncompressibleHyperelasticity::assembleMixedDirichletDeformationDependentNeumann_homogeneousForcing
void assembleMixedDirichletDeformationDependentNeumann_homogeneousForcing(Epetra_Vector &x, Epetra_FECrsMatrix &K, Epetra_FEVector &F)
Definition: nearlyIncompressibleHyperelasticity.cpp:29
neumannInnerSurface_PolyconvexHGO::mu1
double mu1
Definition: neumannInnerSurface_PolyconvexHGO.hpp:18
i
for i
Definition: costFunction.m:38
baseClassFEM::OverlapMap
Epetra_Map * OverlapMap
Definition: baseClassFEM.hpp:17
neumannInnerSurface_PolyconvexHGO::apply_dirichlet_conditions
void apply_dirichlet_conditions(Epetra_FECrsMatrix &K, Epetra_FEVector &F, double &displacement)
Definition: neumannInnerSurface_PolyconvexHGO.hpp:87
neumannInnerSurface_PolyconvexHGO::setup_dirichlet_conditions
void setup_dirichlet_conditions()
Definition: neumannInnerSurface_PolyconvexHGO.hpp:80
laplace::solve_aztec
void solve_aztec(Teuchos::ParameterList &Parameters, Epetra_FECrsMatrix &matrix, Epetra_Vector &lhs, Epetra_FEVector &rhs, int *bc_indx, double *bc_val)
Definition: laplacepp.cpp:81
neumannInnerSurface_PolyconvexHGO::beta4
double beta4
Definition: neumannInnerSurface_PolyconvexHGO.hpp:18
neumannInnerSurface_PolyconvexHGO::setup_clamp
void setup_clamp()
Definition: neumannInnerSurface_PolyconvexHGO.hpp:201