Trilinos based (stochastic) FEM solvers
ceeSBVP.hpp
Go to the documentation of this file.
1 /*
2 Brian Staber (brian.staber@gmail.com)
3 */
4 
5 #ifndef CEESBVP_HPP
6 #define CEESBVP_HPP
7 
8 #include "shinozukapp_layeredcomp_2d.hpp"
9 #include "tensor_calculus.hpp"
10 #include "linearizedElasticity.hpp"
11 
12 class ceeSBVP : public linearizedElasticity
13 {
14 public:
15 
16  Epetra_Vector * solution;
17  Teuchos::ParameterList * Krylov;
18  Teuchos::RCP<shinozuka_layeredcomp_2d> GRF_Generator;
19 
20  Epetra_SerialDenseVector m1, m2, m3, m4, m5;
21  std::vector<int> phase;
22 
23  double _deltaN, _deltaM4, _deltaM5;
24 
25  ceeSBVP(Epetra_Comm & comm, Teuchos::ParameterList & Parameters){
26 
27  Krylov = &Parameters.sublist("Krylov");
28 
29  std::string mesh_file = Teuchos::getParameter<std::string>(Parameters.sublist("Mesh"), "mesh_file");
30  Mesh = new mesh(comm, mesh_file, 1.0);
31  Comm = Mesh->Comm;
32 
33  StandardMap = new Epetra_Map(-1,3*Mesh->n_local_nodes_without_ghosts,&Mesh->local_dof_without_ghosts[0],0,*Comm);
34  OverlapMap = new Epetra_Map(-1,3*Mesh->n_local_nodes,&Mesh->local_dof[0],0,*Comm);
35  ImportToOverlapMap = new Epetra_Import(*OverlapMap,*StandardMap);
36  create_FECrsGraph();
37 
38  setup_dirichlet_conditions();
39  for (unsigned int e=0; e<Mesh->n_cells/32; ++e){
40  for (unsigned int j=0; j<32; ++j){
41  phase.push_back(j);
42  }
43  }
44 
45  int order = Teuchos::getParameter<int>(Parameters.sublist("Shinozuka"), "order");
46  double L1 = Teuchos::getParameter<double>(Parameters.sublist("Shinozuka"), "lx");
47  double L2 = Teuchos::getParameter<double>(Parameters.sublist("Shinozuka"), "ly");
48  _deltaN = Teuchos::getParameter<double>(Parameters.sublist("Shinozuka"), "deltaN");
49  _deltaM4 = Teuchos::getParameter<double>(Parameters.sublist("Shinozuka"), "deltaM4");
50  _deltaM5 = Teuchos::getParameter<double>(Parameters.sublist("Shinozuka"), "deltaM5");
51 
52  GRF_Generator = Teuchos::rcp(new shinozuka_layeredcomp_2d(order,L1,L2));
53 
54  solution = new Epetra_Vector(*StandardMap);
55  }
56 
57  ~ceeSBVP(){
58  }
59 
60  Epetra_SerialDenseVector get_neumannBc(Epetra_SerialDenseMatrix & matrix_X, Epetra_SerialDenseMatrix & xg, unsigned int & gp){
61  std::cout << "Not using this method in this application.\n";
62  Epetra_SerialDenseVector f(3);
63  return f;
64  }
65  Epetra_SerialDenseVector get_forcing(double & x1, double & x2, double & x3, unsigned int & e_lid, unsigned int & gp){
66  std::cout << "Not using this method in this application.\n";
67  Epetra_SerialDenseVector f(3);
68  return f;
69  }
70 
71  void solveOneRealization(double & bcDisp, int * seeds){
72 
73  Epetra_SerialDenseVector w1_shino(Mesh->n_local_cells*Mesh->n_gauss_cells);
74  Epetra_SerialDenseVector w2_shino(Mesh->n_local_cells*Mesh->n_gauss_cells);
75  Epetra_SerialDenseVector w3_shino(Mesh->n_local_cells*Mesh->n_gauss_cells);
76  Epetra_SerialDenseVector w4_shino(Mesh->n_local_cells*Mesh->n_gauss_cells);
77  Epetra_SerialDenseVector w5_shino(Mesh->n_local_cells*Mesh->n_gauss_cells);
78 
79  m1.Resize(Mesh->n_local_cells*Mesh->n_gauss_cells);
80  m2.Resize(Mesh->n_local_cells*Mesh->n_gauss_cells);
81  m3.Resize(Mesh->n_local_cells*Mesh->n_gauss_cells);
82  m4.Resize(Mesh->n_local_cells*Mesh->n_gauss_cells);
83  m5.Resize(Mesh->n_local_cells*Mesh->n_gauss_cells);
84 
85  GRF_Generator->rng.seed(seeds[0]);
86  GRF_Generator->generator_gauss_points(w1_shino,*Mesh,phase);
87 
88  GRF_Generator->rng.seed(seeds[1]);
89  GRF_Generator->generator_gauss_points(w2_shino,*Mesh,phase);
90 
91  GRF_Generator->rng.seed(seeds[2]);
92  GRF_Generator->generator_gauss_points(w3_shino,*Mesh,phase);
93 
94  GRF_Generator->rng.seed(seeds[3]);
95  GRF_Generator->generator_gauss_points(w4_shino,*Mesh,phase);
96 
97  GRF_Generator->rng.seed(seeds[4]);
98  GRF_Generator->generator_gauss_points(w5_shino,*Mesh,phase);
99 
100  double alpha; double beta = 1.0;
101  double Psi1, Psi2;
102 
103  for (unsigned int i=0; i<w1_shino.Length(); ++i){
104  alpha = 3.0/(2.0*_deltaN*_deltaN); beta = 1.0;
105  Psi1 = icdf_gamma(w1_shino(i),alpha,beta);
106  m1(i) = (_deltaN*_deltaN/3.0)*2.0*Psi1;
107 
108  alpha = 3.0/(2.0*_deltaN*_deltaN) - 1.0/2.0;
109  Psi2 = icdf_gamma(w2_shino(i),alpha,beta);
110  m2(i) = (_deltaN*_deltaN/3.0)*( 2.0*Psi2 + w3_shino(i)*w3_shino(i) );
111 
112  m3(i) = (_deltaN*_deltaN/3.0)*std::sqrt(2.0*Psi1)*w3_shino(i);
113 
114  alpha = 1.0/(_deltaM4*_deltaM4); beta = 1.0*_deltaM4*_deltaM4;
115  m4(i) = icdf_gamma(w4_shino(i),alpha,beta);
116  alpha = 1.0/(_deltaM5*_deltaM5); beta = 1.0*_deltaM5*_deltaM5;
117  m5(i) = icdf_gamma(w5_shino(i),alpha,beta);
118  }
119 
120  Epetra_FECrsMatrix linearOperator(Copy,*FEGraph);
121  Epetra_FEVector rhs(*StandardMap);
122  Epetra_Vector lhs(*StandardMap);
123 
124  rhs.PutScalar(0.0);
125  lhs.PutScalar(0.0);
126 
127  assemblePureDirichlet_homogeneousForcing(linearOperator);
128  apply_dirichlet_conditions(linearOperator,rhs,bcDisp);
129 
130  Epetra_LinearProblem problem;
131  AztecOO solver;
132 
133  problem.SetOperator(&linearOperator);
134  problem.SetLHS(&lhs);
135  problem.SetRHS(&rhs);
136  solver.SetProblem(problem);
137  solver.SetParameters(*Krylov);
138 
139  solver.Iterate(2000,1e-6);
140 
141  *solution = lhs;
142  }
143 
144  double icdf_gamma(double & w, double & alpha, double & beta){
145  double erfx = boost::math::erf<double>(w/std::sqrt(2.0));
146  double y = (1.0/2.0)*(1.0 + erfx);
147  double yinv = boost::math::gamma_p_inv<double,double>(alpha,y);
148  double z = yinv*beta;
149  return z;
150  }
151 
152  void setup_dirichlet_conditions(){
153  n_bc_dof = 0;
154  double y;
155  unsigned int node;
156  for (unsigned int i=0; i<Mesh->n_local_nodes_without_ghosts; ++i){
157  node = Mesh->local_nodes[i];
158  y = Mesh->nodes_coord[3*node+1];
159  if(y==0.0){
160  n_bc_dof+=3;
161  }
162  if(y==25.0){
163  n_bc_dof+=3;
164  }
165  }
166 
167  int indbc = 0;
168  dof_on_boundary = new int [n_bc_dof];
169  for (unsigned int inode=0; inode<Mesh->n_local_nodes_without_ghosts; ++inode){
170  node = Mesh->local_nodes[inode];
171  y = Mesh->nodes_coord[3*node+1];
172  if (y==0.0){
173  dof_on_boundary[indbc+0] = 3*inode+0;
174  dof_on_boundary[indbc+1] = 3*inode+1;
175  dof_on_boundary[indbc+2] = 3*inode+2;
176  indbc+=3;
177  }
178  if (y==25.0){
179  dof_on_boundary[indbc+0] = 3*inode+0;
180  dof_on_boundary[indbc+1] = 3*inode+1;
181  dof_on_boundary[indbc+2] = 3*inode+2;
182  indbc+=3;
183  }
184  }
185  }
186 
187  void apply_dirichlet_conditions(Epetra_FECrsMatrix & K, Epetra_FEVector & F, double & displacement){
188  Epetra_MultiVector v(*StandardMap,true);
189  v.PutScalar(0.0);
190 
191  int node;
192  int dof = 1;
193  double y;
194  for (unsigned int inode=0; inode<Mesh->n_local_nodes_without_ghosts; ++inode){
195  node = Mesh->local_nodes[inode];
196  y = Mesh->nodes_coord[3*node+1];
197  if (y==25.0){
198  v[0][StandardMap->LID(3*node+1)] = displacement;
199  }
200  }
201 
202  Epetra_MultiVector rhs_dir(*StandardMap,true);
203  K.Apply(v,rhs_dir);
204  F.Update(-1.0,rhs_dir,1.0);
205 
206  for (unsigned int inode=0; inode<Mesh->n_local_nodes_without_ghosts; ++inode){
207  node = Mesh->local_nodes[inode];
208  y = Mesh->nodes_coord[3*node+1];
209  if (y==0.0){
210  F[0][StandardMap->LID(3*node+0)] = 0.0;
211  F[0][StandardMap->LID(3*node+1)] = 0.0;
212  F[0][StandardMap->LID(3*node+2)] = 0.0;
213  }
214  if (y==25.0){
215  F[0][StandardMap->LID(3*node+0)] = 0.0;
216  F[0][StandardMap->LID(3*node+1)] = displacement;
217  F[0][StandardMap->LID(3*node+2)] = 0.0;
218  }
219  }
220  ML_Epetra::Apply_OAZToMatrix(dof_on_boundary,n_bc_dof,K);
221  }
222 
223  void get_elasticity_tensor(unsigned int & e_lid, unsigned int & gp, Epetra_SerialDenseMatrix & tangent_matrix){
224 
225  int e_gid = Mesh->local_cells[e_lid];
226  int n_gauss_cells = Mesh->n_gauss_cells;
227 
228  double epsilon = 1.0e-6;
229  Epetra_SerialDenseMatrix M(6,6);
230  double M1 = m1(e_lid*n_gauss_cells+gp) + epsilon;
231  double M2 = m2(e_lid*n_gauss_cells+gp) + epsilon;
232  double M3 = m3(e_lid*n_gauss_cells+gp);
233  double M4 = m4(e_lid*n_gauss_cells+gp) + epsilon;
234  double M5 = m5(e_lid*n_gauss_cells+gp) + epsilon;
235 
236  transverse_isotropic_matrix(M,M1,M2,M3,M4,M5);
237 
238  double c1 = 144.8969*1.0e3;
239  double c2 = 14.2500*1.0e3;
240  double c3 = 5.8442*1.0e3;
241  double c4 = 7.5462*1.0e3;
242  double c5 = 12.5580*1.0e3;
243 
244  Epetra_SerialDenseMatrix sqrtmCmoy(6,6);
245  double constant = std::sqrt(c1*c1-2.0*c1*c2+c2*c2+4.0*c3*c3);
246 
247  double d1 = (std::sqrt(2.0)*(c1-c2+constant)*std::sqrt(c1+c2+constant))/(4.0*constant) + (std::sqrt(2.0)*(c2-c1+constant)*std::sqrt(c1+c2-constant))/(4.0*constant);
248  double d2 = (std::sqrt(2.0)*(c2-c1+constant)*std::sqrt(c1+c2+constant))/(4.0*constant) + (std::sqrt(2.0)*(c1-c2+constant)*std::sqrt(c1+c2-constant))/(4.0*constant);
249  double d3 = std::sqrt(2.0)*c3*( std::sqrt(c1+c2+constant) - std::sqrt(c1+c2-constant) )/( 2.0*constant );
250  double d4 = std::sqrt(c4);
251  double d5 = std::sqrt(c5);
252 
253  transverse_isotropic_matrix(sqrtmCmoy,d1,d2,d3,d4,d5);
254 
255  Epetra_SerialDenseMatrix AtimesB(6,6);
256 
257  AtimesB.Multiply('N','N',1.0,M,sqrtmCmoy,0.0);
258  tangent_matrix.Multiply('N','N',1.0,sqrtmCmoy,AtimesB,0.0);
259 
260  if(phase[e_gid] % 2){
261  tangent_matrix(0,5) = -tangent_matrix(0,5);
262  tangent_matrix(5,0) = -tangent_matrix(5,0);
263  tangent_matrix(1,5) = -tangent_matrix(1,5);
264  tangent_matrix(5,1) = -tangent_matrix(5,1);
265  tangent_matrix(2,5) = -tangent_matrix(2,5);
266  tangent_matrix(5,2) = -tangent_matrix(5,2);
267  tangent_matrix(3,4) = -tangent_matrix(3,4);
268  tangent_matrix(4,3) = -tangent_matrix(4,3);
269  }
270  tangent_matrix.Scale(1.0/(1.0+epsilon));
271  }
272 
273  void recover_cauchy_stress(std::string & filename, int * seeds){
274 
275  Epetra_SerialDenseVector w1_shino(Mesh->n_local_cells);
276  Epetra_SerialDenseVector w2_shino(Mesh->n_local_cells);
277  Epetra_SerialDenseVector w3_shino(Mesh->n_local_cells);
278  Epetra_SerialDenseVector w4_shino(Mesh->n_local_cells);
279  Epetra_SerialDenseVector w5_shino(Mesh->n_local_cells);
280 
281  m1.Resize(Mesh->n_local_cells);
282  m2.Resize(Mesh->n_local_cells);
283  m3.Resize(Mesh->n_local_cells);
284  m4.Resize(Mesh->n_local_cells);
285  m5.Resize(Mesh->n_local_cells);
286 
287  double xi = 0.0; double eta = 0.0; double zeta = 0.0;
288 
289  GRF_Generator->rng.seed(seeds[0]);
290  GRF_Generator->generator_one_gauss_point(w1_shino,*Mesh,phase,xi,eta,zeta);
291 
292  GRF_Generator->rng.seed(seeds[1]);
293  GRF_Generator->generator_one_gauss_point(w2_shino,*Mesh,phase,xi,eta,zeta);
294 
295  GRF_Generator->rng.seed(seeds[2]);
296  GRF_Generator->generator_one_gauss_point(w3_shino,*Mesh,phase,xi,eta,zeta);
297 
298  GRF_Generator->rng.seed(seeds[3]);
299  GRF_Generator->generator_one_gauss_point(w4_shino,*Mesh,phase,xi,eta,zeta);
300 
301  GRF_Generator->rng.seed(seeds[4]);
302  GRF_Generator->generator_one_gauss_point(w5_shino,*Mesh,phase,xi,eta,zeta);
303 
304  double alpha; double beta = 1.0;
305  double Psi1, Psi2;
306 
307  for (unsigned int i=0; i<w1_shino.Length(); ++i){
308  alpha = 3.0/(2.0*_deltaN*_deltaN); beta = 1.0;
309  Psi1 = icdf_gamma(w1_shino(i),alpha,beta);
310  m1(i) = (_deltaN*_deltaN/3.0)*2.0*Psi1;
311 
312  alpha = 3.0/(2.0*_deltaN*_deltaN) - 1.0/2.0;
313  Psi2 = icdf_gamma(w2_shino(i),alpha,beta);
314  m2(i) = (_deltaN*_deltaN/3.0)*( 2.0*Psi2 + w3_shino(i)*w3_shino(i) );
315 
316  m3(i) = (_deltaN*_deltaN/3.0)*std::sqrt(2.0*Psi1)*w3_shino(i);
317 
318  alpha = 1.0/(_deltaM4*_deltaM4); beta = 1.0*_deltaM5*_deltaM5;
319  m4(i) = icdf_gamma(w4_shino(i),alpha,beta);
320  alpha = 1.0/(_deltaM4*_deltaM4); beta = 1.0*_deltaM5*_deltaM5;
321  m5(i) = icdf_gamma(w5_shino(i),alpha,beta);
322  }
323  compute_center_cauchy_stress(*solution, filename, false, true);
324  }
325 
326  void get_elasticity_tensor_for_recovery(unsigned int & e_lid, Epetra_SerialDenseMatrix & tangent_matrix){
327 
328  int e_gid = Mesh->local_cells[e_lid];
329 
330  double epsilon = 1.0e-6;
331  Epetra_SerialDenseMatrix M(6,6);
332  double M1 = m1(e_lid) + epsilon;
333  double M2 = m2(e_lid) + epsilon;
334  double M3 = m3(e_lid);
335  double M4 = m4(e_lid) + epsilon;
336  double M5 = m5(e_lid) + epsilon;
337 
338  transverse_isotropic_matrix(M,M1,M2,M3,M4,M5);
339 
340  double c1 = 144.8969*1.0e3;
341  double c2 = 14.2500*1.0e3;
342  double c3 = 5.8442*1.0e3;
343  double c4 = 7.5462*1.0e3;
344  double c5 = 12.5580*1.0e3;
345 
346  Epetra_SerialDenseMatrix sqrtmCmoy(6,6);
347  double constant = std::sqrt(c1*c1-2.0*c1*c2+c2*c2+4.0*c3*c3);
348 
349  double d1 = (std::sqrt(2.0)*(c1-c2+constant)*std::sqrt(c1+c2+constant))/(4.0*constant) + (std::sqrt(2.0)*(c2-c1+constant)*std::sqrt(c1+c2-constant))/(4.0*constant);
350  double d2 = (std::sqrt(2.0)*(c2-c1+constant)*std::sqrt(c1+c2+constant))/(4.0*constant) + (std::sqrt(2.0)*(c1-c2+constant)*std::sqrt(c1+c2-constant))/(4.0*constant);
351  double d3 = std::sqrt(2.0)*c3*( std::sqrt(c1+c2+constant) - std::sqrt(c1+c2-constant) )/( 2.0*constant );
352  double d4 = std::sqrt(c4);
353  double d5 = std::sqrt(c5);
354 
355  transverse_isotropic_matrix(sqrtmCmoy,d1,d2,d3,d4,d5);
356 
357  Epetra_SerialDenseMatrix AtimesB(6,6);
358 
359  AtimesB.Multiply('N','N',1.0,M,sqrtmCmoy,0.0);
360  tangent_matrix.Multiply('N','N',1.0,sqrtmCmoy,AtimesB,0.0);
361 
362  if(phase[e_gid] % 2){
363  tangent_matrix(0,5) = -tangent_matrix(0,5);
364  tangent_matrix(5,0) = -tangent_matrix(5,0);
365  tangent_matrix(1,5) = -tangent_matrix(1,5);
366  tangent_matrix(5,1) = -tangent_matrix(5,1);
367  tangent_matrix(2,5) = -tangent_matrix(2,5);
368  tangent_matrix(5,2) = -tangent_matrix(5,2);
369  tangent_matrix(3,4) = -tangent_matrix(3,4);
370  tangent_matrix(4,3) = -tangent_matrix(4,3);
371  }
372  tangent_matrix.Scale(1.0/(1.0+epsilon));
373  }
374 
375  void transverse_isotropic_matrix(Epetra_SerialDenseMatrix & C, double & c1, double & c2, double & c3, double & c4, double & c5){
376 
377  C(0,0) = c1/16.0 + (9.0*c2)/32.0 + (9.0*c4)/32.0 + (3.0*c5)/8.0 + (3.0*std::sqrt(2.0)*c3)/16.0;
378  C(0,1) = (3.0*c1)/16.0 + (3.0*c2)/32.0 + (3.0*c4)/32.0 - (3.0*c5)/8.0 + (5.0*std::sqrt(2.0)*c3)/16.0;
379  C(0,2) = (3.0*c2)/8.0 - (3.0*c4)/8.0 + (std::sqrt(2.0)*c3)/8.0;
380  C(0,3) = 0.0;
381  C(0,4) = 0.0;
382  C(0,5) = std::sqrt(6)*(c1/16.0 - (3.0*c2)/32.0 - (3.0*c4)/32.0 + c5/8.0 + (std::sqrt(2.0)*c3)/16.0);
383 
384  C(1,0) = (3.0*c1)/16.0 + (3.0*c2)/32.0 + (3.0*c4)/32.0 - (3.0*c5)/8.0 + (5*std::sqrt(2.0)*c3)/16.0;
385  C(1,1) = (9.0*c1)/16.0 + c2/32.0 + c4/32.0 + (3.0*c5)/8.0 + (3.0*std::sqrt(2.0)*c3)/16.0;
386  C(1,2) = c2/8.0 - c4/8.0 + (3.0*std::sqrt(2.0)*c3)/8.0;
387  C(1,3) = 0.0;
388  C(1,4) = 0.0;
389  C(1,5) = -std::sqrt(6.0)*(c2/32.0 - (3.0*c1)/16.0 + c4/32.0 + c5/8.0 + (std::sqrt(2.0)*c3)/16.0);
390 
391  C(2,0) = (3.0*c2)/8.0 - (3.0*c4)/8.0 + std::sqrt(2.0)*c3/8.0;
392  C(2,1) = c2/8.0 - c4/8.0 + (3.0*std::sqrt(2.0)*c3)/8.0;
393  C(2,2) = c2/2.0 + c4/2.0;
394  C(2,3) = 0.0;
395  C(2,4) = 0.0;
396  C(2,5) = (std::sqrt(3.0)*(2.0*c3 - std::sqrt(2.0)*c2 + std::sqrt(2.0)*c4))/8.0;
397 
398  C(3,0) = 0.0;
399  C(3,1) = 0.0;
400  C(3,2) = 0.0;
401  C(3,3) = c4/4.0 + (3.0*c5)/4.0;
402  C(3,4) = -(std::sqrt(3.0)*(c4 - c5))/4.0;
403  C(3,5) = 0.0;
404 
405  C(4,0) = 0.0;
406  C(4,1) = 0.0;
407  C(4,2) = 0.0;
408  C(4,3) = -(std::sqrt(3.0)*(c4 - c5))/4.0;
409  C(4,4) = (3.0*c4)/4.0 + c5/4.0;
410  C(4,5) = 0.0;
411 
412  C(5,0) = std::sqrt(6.0)*(c1/16.0 - (3.0*c2)/32.0 - (3.0*c4)/32.0 + c5/8.0 + std::sqrt(2.0)*c3/16.0);
413  C(5,1) = -std::sqrt(6.0)*(c2/32.0 - (3.0*c1)/16.0 + c4/32.0 + c5/8.0 + (std::sqrt(2.0)*c3)/16.0);
414  C(5,2) = (std::sqrt(3.0)*(2.0*c3 - std::sqrt(2.0)*c2 + std::sqrt(2.0)*c4))/8.0;
415  C(5,3) = 0.0;
416  C(5,4) = 0.0;
417  C(5,5) = (3.0*c1)/8.0 + (3.0*c2)/16.0 + (3.0*c4)/16.0 + c5/4.0 - (3.0*std::sqrt(2.0)*c3)/8.0;
418 
419  }
420 
421  int print_solution(std::string fileName){
422 
423  int NumTargetElements = 0;
424  if (Comm->MyPID()==0){
425  NumTargetElements = 3*Mesh->n_nodes;
426  }
427  Epetra_Map MapOnRoot(-1,NumTargetElements,0,*Comm);
428  Epetra_Export ExportOnRoot(*StandardMap,MapOnRoot);
429  Epetra_MultiVector lhs_root(MapOnRoot,true);
430  lhs_root.Export(*solution,ExportOnRoot,Insert);
431 
432  int error = EpetraExt::MultiVectorToMatrixMarketFile(fileName.c_str(),lhs_root,0,0,false);
433 
434  return error;
435 
436  }
437 };
438 
439 
440 #endif
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
ceeSBVP::Krylov
Teuchos::ParameterList * Krylov
Definition: ceeSBVP.hpp:17
mesh::nodes_coord
std::vector< double > nodes_coord
Definition: meshpp.hpp:77
ceeSBVP::ceeSBVP
ceeSBVP(Epetra_Comm &comm, Teuchos::ParameterList &Parameters)
Definition: ceeSBVP.hpp:25
ceeSBVP::_deltaM4
double _deltaM4
Definition: ceeSBVP.hpp:23
ceeSBVP::m4
Epetra_SerialDenseVector m4
Definition: ceeSBVP.hpp:20
ceeSBVP::GRF_Generator
Teuchos::RCP< shinozuka_layeredcomp_2d > GRF_Generator
Definition: ceeSBVP.hpp:18
ceeSBVP::solveOneRealization
void solveOneRealization(double &bcDisp, int *seeds)
Definition: ceeSBVP.hpp:71
mesh::n_local_cells
int n_local_cells
Definition: meshpp.hpp:94
ceeSBVP::solution
Epetra_Vector * solution
Definition: ceeSBVP.hpp:16
ceeSBVP::recover_cauchy_stress
void recover_cauchy_stress(std::string &filename, int *seeds)
Definition: ceeSBVP.hpp:273
ceeSBVP::phase
std::vector< int > phase
Definition: ceeSBVP.hpp:21
ceeSBVP::m2
Epetra_SerialDenseVector m2
Definition: ceeSBVP.hpp:20
ceeSBVP::~ceeSBVP
~ceeSBVP()
Definition: ceeSBVP.hpp:57
ceeSBVP::apply_dirichlet_conditions
void apply_dirichlet_conditions(Epetra_FECrsMatrix &K, Epetra_FEVector &F, double &displacement)
Definition: ceeSBVP.hpp:187
baseClassFEM::FEGraph
Epetra_FECrsGraph * FEGraph
Definition: baseClassFEM.hpp:20
ceeSBVP::icdf_gamma
double icdf_gamma(double &w, double &alpha, double &beta)
Definition: ceeSBVP.hpp:144
ceeSBVP::print_solution
int print_solution(std::string fileName)
Definition: ceeSBVP.hpp:421
ceeSBVP::get_elasticity_tensor_for_recovery
void get_elasticity_tensor_for_recovery(unsigned int &e_lid, Epetra_SerialDenseMatrix &tangent_matrix)
Definition: ceeSBVP.hpp:326
mesh::local_cells
std::vector< int > local_cells
Definition: meshpp.hpp:81
mesh
clc clear all close all mesh
Definition: run.m:5
mesh::n_local_nodes_without_ghosts
int n_local_nodes_without_ghosts
Definition: meshpp.hpp:92
filename
filename
Definition: costFunction.m:44
linearizedElasticity::compute_center_cauchy_stress
void compute_center_cauchy_stress(Epetra_Vector &x, std::string &filename, bool printCauchy, bool printVM)
Definition: linearizedElasticity.cpp:433
linearizedElasticity::assemblePureDirichlet_homogeneousForcing
void assemblePureDirichlet_homogeneousForcing(Epetra_FECrsMatrix &K)
Definition: linearizedElasticity.cpp:56
baseClassFEM::ImportToOverlapMap
Epetra_Import * ImportToOverlapMap
Definition: baseClassFEM.hpp:19
baseClassFEM::StandardMap
Epetra_Map * StandardMap
Definition: baseClassFEM.hpp:18
ceeSBVP::m5
Epetra_SerialDenseVector m5
Definition: ceeSBVP.hpp:20
mesh::n_gauss_cells
unsigned int n_gauss_cells
Definition: meshpp.hpp:98
ceeSBVP::transverse_isotropic_matrix
void transverse_isotropic_matrix(Epetra_SerialDenseMatrix &C, double &c1, double &c2, double &c3, double &c4, double &c5)
Definition: ceeSBVP.hpp:375
mesh::n_nodes
int n_nodes
Definition: meshpp.hpp:87
beta
modelParameters beta
Definition: run_station15.m:10
i
for i
Definition: costFunction.m:38
baseClassFEM::OverlapMap
Epetra_Map * OverlapMap
Definition: baseClassFEM.hpp:17
ceeSBVP::get_elasticity_tensor
void get_elasticity_tensor(unsigned int &e_lid, unsigned int &gp, Epetra_SerialDenseMatrix &tangent_matrix)
Definition: ceeSBVP.hpp:223
e
e
Definition: run.m:10
ceeSBVP::setup_dirichlet_conditions
void setup_dirichlet_conditions()
Definition: ceeSBVP.hpp:152
shinozukapp_layeredcomp_2d.hpp
ceeSBVP::_deltaN
double _deltaN
Definition: ceeSBVP.hpp:23
ceeSBVP::get_forcing
Epetra_SerialDenseVector get_forcing(double &x1, double &x2, double &x3, unsigned int &e_lid, unsigned int &gp)
Definition: ceeSBVP.hpp:65
shinozuka_layeredcomp_2d
Definition: shinozukapp_layeredcomp_2d.hpp:16
ceeSBVP::get_neumannBc
Epetra_SerialDenseVector get_neumannBc(Epetra_SerialDenseMatrix &matrix_X, Epetra_SerialDenseMatrix &xg, unsigned int &gp)
Definition: ceeSBVP.hpp:60
eta
output eta
Definition: costFunction.m:33
ceeSBVP::m1
Epetra_SerialDenseVector m1
Definition: ceeSBVP.hpp:20
ceeSBVP::_deltaM5
double _deltaM5
Definition: ceeSBVP.hpp:23
ceeSBVP::m3
Epetra_SerialDenseVector m3
Definition: ceeSBVP.hpp:20