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