My Project
elasticity_upscale_impl.hpp
Go to the documentation of this file.
1 //==============================================================================
11 //==============================================================================
12 #ifndef OPM_ELASTICITY_UPSCALE_IMPL_HPP
13 #define OPM_ELASTICITY_UPSCALE_IMPL_HPP
14 
15 #include <iostream>
16 
17 #ifdef HAVE_OPENMP
18 #include <omp.h>
19 #endif
20 
21 #include <opm/parser/eclipse/Deck/DeckKeyword.hpp>
22 
23 namespace Opm {
24 namespace Elasticity {
25 
26 #undef IMPL_FUNC
27 #define IMPL_FUNC(A,B) template<class GridType, class PC> \
28  A ElasticityUpscale<GridType, PC>::B
29 
30 IMPL_FUNC(std::vector<BoundaryGrid::Vertex>,
31  extractFace(Direction dir, ctype coord))
32 {
33  std::vector<BoundaryGrid::Vertex> result;
34  const LeafVertexIterator itend = gv.leafGridView().template end<dim>();
35 
36  // make a mapper for codim dim entities in the leaf grid
37  using LeafGridView = Dune::GridView<Dune::DefaultLeafGridViewTraits<GridType>>;
38  Dune::MultipleCodimMultipleGeomTypeMapper<LeafGridView> mapper(gv.leafGridView(), Dune::mcmgVertexLayout());
39  // iterate over vertices and find slaves
40  LeafVertexIterator start = gv.leafGridView().template begin<dim>();
41  for (LeafVertexIterator it = start; it != itend; ++it) {
42  if (isOnPlane(dir,it->geometry().corner(0),coord)) {
43  BoundaryGrid::Vertex v;
44  v.i = mapper.index(*it);
45  BoundaryGrid::extract(v.c,it->geometry().corner(0),log2(float(dir)));
46  result.push_back(v);
47  }
48  }
49 
50  return result;
51 }
52 
53 
54 IMPL_FUNC(BoundaryGrid, extractMasterFace(Direction dir,
55  ctype coord,
56  SIDE side,
57  bool dc))
58 {
59  static const int V1[3][4] = {{0,2,4,6},
60  {0,1,4,5},
61  {0,1,2,3}};
62  static const int V2[3][4] = {{1,3,5,7},
63  {2,3,6,7},
64  {4,5,6,7}};
65  const LeafIndexSet& set = gv.leafGridView().indexSet();
66 
67  int c = 0;
68  int i = log2(float(dir));
69  BoundaryGrid result;
70  // we first group nodes into this map through the coordinate of lower left
71  // vertex. we then split this up into pillars for easy processing later
72  std::map<double, std::vector<BoundaryGrid::Quad> > nodeMap;
73  for (LeafIterator cell = gv.leafGridView().template begin<0>();
74  cell != gv.leafGridView().template end<0>(); ++cell, ++c) {
75  std::vector<BoundaryGrid::Vertex> verts;
76  int idx=0;
77  if (side == LEFT)
78  idx = set.subIndex(*cell,V1[i][0],dim);
79  else if (side == RIGHT)
80  idx = set.subIndex(*cell,V2[i][0],dim);
81  Dune::FieldVector<double, 3> pos = gv.vertexPosition(idx);
82  if (isOnPlane(dir,pos,coord)) {
83  for (int j=0;j<4;++j) {
84  if (side == LEFT)
85  idx = set.subIndex(*cell,V1[i][j],dim);
86  if (side == RIGHT)
87  idx = set.subIndex(*cell,V2[i][j],dim);
88  pos = gv.vertexPosition(idx);
89  if (!isOnPlane(dir,pos,coord))
90  continue;
91  BoundaryGrid::Vertex v;
92  BoundaryGrid::extract(v,pos,i);
93  v.i = idx;
94  verts.push_back(v);
95  }
96  }
97  if (verts.size() == 4) {
98  BoundaryGrid::Quad q;
99  q.v[0] = minXminY(verts);
100  q.v[1] = maxXminY(verts);
101  if (dc) {
102  q.v[2] = minXmaxY(verts);
103  q.v[3] = maxXmaxY(verts);
104  } else {
105  q.v[2] = maxXmaxY(verts);
106  q.v[3] = minXmaxY(verts);
107  }
108  std::map<double, std::vector<BoundaryGrid::Quad> >::iterator it;
109  for (it = nodeMap.begin(); it != nodeMap.end(); ++it) {
110  if (fabs(it->first-q.v[0].c[0]) < 1.e-7) {
111  it->second.push_back(q);
112  break;
113  }
114  }
115  if (it == nodeMap.end())
116  nodeMap[q.v[0].c[0]].push_back(q);
117 
118  result.add(q);
119  }
120  }
121 
122  int p=0;
123  std::map<double, std::vector<BoundaryGrid::Quad> >::const_iterator it;
124  for (it = nodeMap.begin(); it != nodeMap.end(); ++it, ++p) {
125  for (size_t ii=0;ii<it->second.size();++ii)
126  result.addToColumn(p,it->second[ii]);
127  }
128 
129  return result;
130 }
131 
132 IMPL_FUNC(void, determineSideFaces(const double* min, const double* max))
133 {
134  master.push_back(extractMasterFace(X,min[0]));
135  master.push_back(extractMasterFace(Y,min[1]));
136  master.push_back(extractMasterFace(Z,min[2]));
137 
138  slave.push_back(extractFace(X,max[0]));
139  slave.push_back(extractFace(Y,max[1]));
140  slave.push_back(extractFace(Z,max[2]));
141 }
142 
143 IMPL_FUNC(void, findBoundaries(double* min, double* max))
144 {
145  max[0] = max[1] = max[2] = -1e5;
146  min[0] = min[1] = min[2] = 1e5;
147  const LeafVertexIterator itend = gv.leafGridView().template end<dim>();
148 
149  // iterate over vertices and find slaves
150  LeafVertexIterator start = gv.leafGridView().template begin<dim>();
151  for (LeafVertexIterator it = start; it != itend; ++it) {
152  for (int i=0;i<3;++i) {
153  min[i] = std::min(min[i],it->geometry().corner(0)[i]);
154  max[i] = std::max(max[i],it->geometry().corner(0)[i]);
155  }
156  }
157 }
158 
159 IMPL_FUNC(void, addMPC(Direction dir, int slavenode,
160  const BoundaryGrid::Vertex& m))
161 {
162  MPC* mpc = new MPC(slavenode,log2(float(dir))+1);
163  if (m.i > -1) { // we matched a node exactly
164  mpc->addMaster(m.i,log2(float(dir))+1,1.f);
165  } else {
166  std::vector<double> N = m.q->evalBasis(m.c[0],m.c[1]);
167  for (int i=0;i<4;++i)
168  mpc->addMaster(m.q->v[i].i,log2(float(dir))+1,N[i]);
169  }
170  A.addMPC(mpc);
171 }
172 
173 IMPL_FUNC(void, periodicPlane(Direction /*plane */,
174  Direction /* dir */,
175  const std::vector<BoundaryGrid::Vertex>& slavepointgrid,
176  const BoundaryGrid& mastergrid))
177 {
178  for (size_t i=0;i<slavepointgrid.size();++i) {
179  BoundaryGrid::Vertex coord;
180  if (mastergrid.find(coord,slavepointgrid[i])) {
181  addMPC(X,slavepointgrid[i].i,coord);
182  addMPC(Y,slavepointgrid[i].i,coord);
183  addMPC(Z,slavepointgrid[i].i,coord);
184  }
185  }
186 }
187 
193 static std::vector< std::vector<int> > renumber(const BoundaryGrid& b,
194  int n1, int n2,
195  int& totalDOFs)
196 {
197  std::vector<std::vector<int> > nodes;
198  nodes.resize(b.size());
199  // loop over elements
200  totalDOFs = 0;
201 
202  // fix lower left multiplicator.
203  // will be "transfered" to all corners through periodic conditions
204  nodes[0].push_back(-1);
205  for (size_t e=0; e < b.size(); ++e) {
206  // first direction major ordered nodes within each element
207  for (int i2=0; i2 < n2; ++i2) {
208  if (e != 0)
209  nodes[e].push_back(nodes[e-1][i2*n1+n1-1]);
210 
211  int start = (e==0 && i2 != 0)?0:1;
212 
213  // slave the buggers
214  if (i2 == n2-1 && n2 > 2) {
215  for (int i1=(e==0?0:1); i1 < n1; ++i1) {
216  nodes[e].push_back(nodes[e][i1]);
217  }
218  } else {
219  for (int i1=start; i1 < n1; ++i1) {
220  if (e == b.size()-1)
221  nodes[e].push_back(nodes[0][i2*n1]);
222  else
223  nodes[e].push_back(totalDOFs++);
224  }
225  }
226  }
227  }
228 
229  return nodes;
230 }
231 
232 IMPL_FUNC(int, addBBlockMortar(const BoundaryGrid& b1,
233  const BoundaryGrid& interface,
234  int /* dir */, int n1, int n2,
235  int colofs))
236 {
237  // renumber the linear grid to the real multiplier grid
238  int totalEqns;
239  std::vector<std::vector<int> > lnodes = renumber(interface, n1+1,
240  n2+1, totalEqns);
241  if (Bpatt.empty())
242  Bpatt.resize(A.getEqns());
243 
244  // process pillar by pillar
245  for (size_t p=0;p<interface.size();++p) {
246  for (size_t q=0;q<b1.colSize(p);++q) {
247  for (size_t i=0;i<4;++i) {
248  for (size_t d=0;d<3;++d) {
249  const MPC* mpc = A.getMPC(b1.getQuad(p,q).v[i].i,d);
250  if (mpc) {
251  for (size_t n=0;n<mpc->getNoMaster();++n) {
252  int dof = A.getEquationForDof(mpc->getMaster(n).node,d);
253  if (dof > -1) {
254  for (size_t j=0;j<lnodes[p].size();++j) {
255  int indexj = 3*lnodes[p][j]+d;
256  if (indexj > -1)
257  Bpatt[dof].insert(indexj+colofs);
258  }
259  }
260  }
261  } else {
262  int dof = A.getEquationForDof(b1.getQuad(p,q).v[i].i,d);
263  if (dof > -1) {
264  for (size_t j=0;j<lnodes[p].size();++j) {
265  int indexj = 3*lnodes[p][j]+d;
266  if (indexj > -1)
267  Bpatt[dof].insert(indexj+colofs);
268  }
269  }
270  }
271  }
272  }
273  }
274  }
275 
276  return 3*totalEqns;
277 }
278 
279 IMPL_FUNC(void, assembleBBlockMortar(const BoundaryGrid& b1,
280  const BoundaryGrid& interface,
281  int dir, int n1,
282  int n2, int colofs,
283  double alpha))
284 {
285  // get a set of P1 shape functions for the displacements
286  P1ShapeFunctionSet<ctype,ctype,2> ubasis =
288 
289  // get a set of PN shape functions for the multipliers
290  PNShapeFunctionSet<2> lbasis(n1+1, n2+1);
291 
292  // renumber the linear grid to the real multiplier grid
293  int totalEqns;
294  std::vector<std::vector<int> > lnodes = renumber(interface, n1+1,
295  n2+1, totalEqns);
296  // get a reference element
297  auto gt = Dune::GeometryTypes::cube(2);
298  // get a quadrature rule
299  int quadorder = std::max((1.0+n1+0.5)/2.0,(1.0+n2+0.5)/2.0);
300  quadorder = std::max(quadorder, 2);
301  const Dune::QuadratureRule<ctype,2>& rule =
302  Dune::QuadratureRules<ctype,2>::rule(gt,quadorder);
303 
304  // do the assembly loop
305  typename Dune::QuadratureRule<ctype,2>::const_iterator r;
306  Dune::DynamicMatrix<ctype> lE(ubasis.n,(n1+1)*(n2+1),0.0);
307  LoggerHelper help(interface.size(), 5, 1000);
308  for (int g=0;g<5;++g) {
309  for (int p=help.group(g).first;p<help.group(g).second;++p) {
310  const BoundaryGrid::Quad& qi(interface[p]);
311  HexGeometry<2,2,GridType> lg(qi);
312  for (size_t q=0;q<b1.colSize(p);++q) {
313  const BoundaryGrid::Quad& qu = b1.getQuad(p,q);
314  HexGeometry<2,2,GridType> hex(qu,gv,dir);
315  lE = 0;
316  for (r = rule.begin(); r != rule.end();++r) {
317  ctype detJ = hex.integrationElement(r->position());
318  if (detJ < 0)
319  assert(0);
320 
321  typename HexGeometry<2,2,GridType>::LocalCoordinate loc =
322  lg.local(hex.global(r->position()));
323  assert(loc[0] <= 1.0+1.e-4 && loc[0] >= 0.0 && loc[1] <= 1.0+1.e-4 && loc[1] >= 0.0);
324  for (int i=0;i<ubasis.n;++i) {
325  for (int j=0;j<lbasis.size();++j) {
326  lE[i][j] += ubasis[i].evaluateFunction(r->position())*
327  lbasis[j].evaluateFunction(loc)*detJ*r->weight();
328  }
329  }
330  }
331 
332  // and assemble element contributions
333  for (int d=0;d<3;++d) {
334  for (int i=0;i<4;++i) {
335  const MPC* mpc = A.getMPC(qu.v[i].i,d);
336  if (mpc) {
337  for (size_t n=0;n<mpc->getNoMaster();++n) {
338  int indexi = A.getEquationForDof(mpc->getMaster(n).node,d);
339  if (indexi > -1) {
340  for (size_t j=0;j<lnodes[p].size();++j) {
341  int indexj = lnodes[p][j]*3+d;
342  if (indexj > -1)
343  B[indexi][indexj+colofs] += alpha*lE[i][j];
344  }
345  }
346  }
347  } else {
348  int indexi = A.getEquationForDof(qu.v[i].i,d);
349  if (indexi > -1) {
350  for (size_t j=0;j<lnodes[p].size();++j) {
351  int indexj = lnodes[p][j]*3+d;
352  if (indexj > -1)
353  B[indexi][indexj+colofs] += alpha*lE[i][j];
354  }
355  }
356  }
357  }
358  }
359  }
360  }
361  help.log(g, "\t\t\t... still processing ... pillar ");
362  }
363 }
364 
365 IMPL_FUNC(void, fixPoint(Direction dir,
366  GlobalCoordinate coord,
367  const NodeValue& value))
368 {
369  typedef typename GridType::LeafGridView::template Codim<dim>::Iterator VertexLeafIterator;
370  const VertexLeafIterator itend = gv.leafGridView().template end<dim>();
371 
372  // make a mapper for codim 0 entities in the leaf grid
373  using LeafGridView = Dune::GridView<Dune::DefaultLeafGridViewTraits<GridType>>;
374  Dune::MultipleCodimMultipleGeomTypeMapper<LeafGridView> mapper(gv.leafGridView(), Dune::mcmgVertexLayout());
375 
376  // iterate over vertices
377  for (VertexLeafIterator it = gv.leafGridView().template begin<dim>(); it != itend; ++it) {
378  if (isOnPoint(it->geometry().corner(0),coord)) {
379  int indexi = mapper.index(*it);
380  A.updateFixedNode(indexi,std::make_pair(dir,value));
381  }
382  }
383 }
384 
385 IMPL_FUNC(bool, isOnPlane(Direction plane,
386  GlobalCoordinate coord,
387  ctype value))
388 {
389  if (plane < X || plane > Z)
390  return false;
391  int p = log2(float(plane));
392  ctype delta = fabs(value-coord[p]);
393  return delta < tol;
394 }
395 
396 IMPL_FUNC(void, fixLine(Direction dir,
397  ctype x, ctype y,
398  const NodeValue& value))
399 {
400  typedef typename GridType::LeafGridView::template Codim<dim>::Iterator VertexLeafIterator;
401  const VertexLeafIterator itend = gv.leafGridView().template end<dim>();
402 
403  // make a mapper for codim 0 entities in the leaf grid
404  using LeafGridView = Dune::GridView<Dune::DefaultLeafGridViewTraits<GridType>>;
405  Dune::MultipleCodimMultipleGeomTypeMapper<LeafGridView> mapper(gv, Dune::mcmgVertexLayout());
406 
407  // iterate over vertices
408  for (VertexLeafIterator it = gv.leafGridView().template begin<dim>(); it != itend; ++it) {
409  if (isOnLine(dir,it->geometry().corner(0),x,y)) {
410  int indexi = mapper.index(*it);
411  A.updateFixedNode(indexi,std::make_pair(XYZ,value));
412  }
413  }
414 }
415 
416 IMPL_FUNC(bool, isOnLine(Direction dir,
417  GlobalCoordinate coord,
418  ctype x, ctype y))
419 {
420  if (dir < X || dir > Z)
421  return false;
422  int ix = int(log2(float(dir))+1) % 3;
423  int iy = int(log2(float(dir))+2) % 3;
424  ctype delta = x-coord[ix];
425  if (delta > tol || delta < -tol)
426  return false;
427  delta = y-coord[iy];
428  if (delta > tol || delta < -tol)
429  return false;
430 
431  return true;
432 }
433 
434 IMPL_FUNC(bool, isOnPoint(GlobalCoordinate coord,
435  GlobalCoordinate point))
436 {
437  GlobalCoordinate delta = point-coord;
438  return delta.one_norm() < tol;
439 }
440 
441 IMPL_FUNC(void, assemble(int loadcase, bool matrix))
442 {
443  const int comp = 3+(dim-2)*3;
444  static const int bfunc = 4+(dim-2)*4;
445 
446  Dune::FieldVector<ctype,comp> eps0;
447  eps0 = 0;
448  if (loadcase > -1) {
449  eps0[loadcase] = 1;
450  b[loadcase] = 0;
451  }
452  if (matrix)
453  A.getOperator() = 0;
454 
455  for (int i=0;i<2;++i) {
456  if (color[1].size() && matrix)
457  std::cout << "\tprocessing " << (i==0?"red ":"black ") << "elements" << std::endl;
458 #pragma omp parallel for schedule(static)
459  for (size_t j=0;j<color[i].size();++j) {
460  Dune::FieldMatrix<ctype,comp,comp> C;
461  Dune::FieldMatrix<ctype,dim*bfunc,dim*bfunc> K;
462  Dune::FieldMatrix<ctype,dim*bfunc,dim*bfunc>* KP=0;
463  Dune::FieldVector<ctype,dim*bfunc> ES;
464  Dune::FieldVector<ctype,dim*bfunc>* EP=0;
465  if (matrix)
466  KP = &K;
467  if (loadcase > -1)
468  EP = &ES;
469 
470  for (size_t k=0;k<color[i][j].size();++k) {
471  LeafIterator it = gv.leafGridView().template begin<0>();
472  for (int l=0;l<color[i][j][k];++l)
473  ++it;
474  materials[color[i][j][k]]->getConstitutiveMatrix(C);
475  // determine geometry type of the current element and get the matching reference element
476  Dune::GeometryType gt = it->type();
477 
478  Dune::FieldMatrix<ctype,dim*bfunc,dim*bfunc> Aq;
479  K = 0;
480  ES = 0;
481 
482  // get a quadrature rule of order two for the given geometry type
483  const Dune::QuadratureRule<ctype,dim>& rule = Dune::QuadratureRules<ctype,dim>::rule(gt,2);
484  for (typename Dune::QuadratureRule<ctype,dim>::const_iterator r = rule.begin();
485  r != rule.end() ; ++r) {
486  // compute the jacobian inverse transposed to transform the gradients
487  Dune::FieldMatrix<ctype,dim,dim> jacInvTra =
488  it->geometry().jacobianInverseTransposed(r->position());
489 
490  ctype detJ = it->geometry().integrationElement(r->position());
491  if (detJ <= 1.e-5 && verbose) {
492  std::cout << "cell " << color[i][j][k] << " is (close to) degenerated, detJ " << detJ << std::endl;
493  double zdiff=0.0;
494  for (int ii=0;ii<4;++ii)
495  zdiff = std::max(zdiff, it->geometry().corner(ii+4)[2]-it->geometry().corner(ii)[2]);
496  std::cout << " - Consider setting ctol larger than " << zdiff << std::endl;
497  }
498 
499  Dune::FieldMatrix<ctype,comp,dim*bfunc> lB;
500  E.getBmatrix(lB,r->position(),jacInvTra);
501 
502  if (matrix) {
503  E.getStiffnessMatrix(Aq,lB,C,detJ*r->weight());
504  K += Aq;
505  }
506 
507  // load vector
508  if (EP) {
509  Dune::FieldVector<ctype,dim*bfunc> temp;
510  temp = Dune::FMatrixHelp::multTransposed(lB,Dune::FMatrixHelp::mult(C,eps0));
511  temp *= -detJ*r->weight();
512  ES += temp;
513  }
514  }
515  A.addElement(KP,EP,it,(loadcase > -1)?&b[loadcase]:NULL);
516  }
517  }
518  }
519 }
520 
521 IMPL_FUNC(template<int comp> void,
522  averageStress(Dune::FieldVector<ctype,comp>& sigma,
523  const Vector& uarg, int loadcase))
524 {
525  if (loadcase < 0 || loadcase > 5)
526  return;
527 
528  static const int bfunc = 4+(dim-2)*4;
529 
530  const LeafIterator itend = gv.leafGridView().template end<0>();
531 
532  Dune::FieldMatrix<ctype,comp,comp> C;
533  Dune::FieldVector<ctype,comp> eps0;
534  eps0 = 0;
535  eps0[loadcase] = 1;
536  int m=0;
537  sigma = 0;
538  double volume=0;
539  for (LeafIterator it = gv.leafGridView().template begin<0>(); it != itend; ++it) {
540  materials[m++]->getConstitutiveMatrix(C);
541  // determine geometry type of the current element and get the matching reference element
542  Dune::GeometryType gt = it->type();
543 
544  Dune::FieldVector<ctype,bfunc*dim> v;
545  A.extractValues(v,uarg,it);
546 
547  // get a quadrature rule of order two for the given geometry type
548  const Dune::QuadratureRule<ctype,dim>& rule = Dune::QuadratureRules<ctype,dim>::rule(gt,2);
549  for (typename Dune::QuadratureRule<ctype,dim>::const_iterator r = rule.begin();
550  r != rule.end() ; ++r) {
551  // compute the jacobian inverse transposed to transform the gradients
552  Dune::FieldMatrix<ctype,dim,dim> jacInvTra =
553  it->geometry().jacobianInverseTransposed(r->position());
554 
555  ctype detJ = it->geometry().integrationElement(r->position());
556 
557  volume += detJ*r->weight();
558 
559  Dune::FieldMatrix<ctype,comp,dim*bfunc> lB;
560  E.getBmatrix(lB,r->position(),jacInvTra);
561 
562  Dune::FieldVector<ctype,comp> s;
563  E.getStressVector(s,v,eps0,lB,C);
564  s *= detJ*r->weight();
565  sigma += s;
566  }
567  }
568  sigma /= volume;
569  if (Escale > 0)
570  sigma /= Escale/Emin;
571 }
572 
573 IMPL_FUNC(void, loadMaterialsFromGrid(const std::string& file))
574 {
575  typedef std::map<std::pair<double,double>,
576  std::shared_ptr<Material> > MaterialMap;
577  MaterialMap cache;
578  std::vector<double> Emod;
579  std::vector<double> Poiss;
580  std::vector<int> satnum;
581  std::vector<double> rho;
582  upscaledRho = -1;
583 
584  if (file == "uniform") {
585  int cells = gv.size(0);
586  Emod.insert(Emod.begin(),cells,100.f);
587  Poiss.insert(Poiss.begin(),cells,0.38f);
588  } else {
589  Opm::Parser parser;
590  auto deck = parser.parseFile(file);
591  if (deck.hasKeyword("YOUNGMOD") && deck.hasKeyword("POISSONMOD")) {
592  Emod = deck.getKeyword("YOUNGMOD").getRawDoubleData();
593  Poiss = deck.getKeyword("POISSONMOD").getRawDoubleData();
594  std::vector<double>::const_iterator it = std::min_element(Poiss.begin(), Poiss.end());
595  if (*it < 0) {
596  std::cerr << "Auxetic material specified for cell " << it-Poiss.begin() << std::endl
597  << "Emod: "<< Emod[it-Poiss.begin()] << " Poisson's ratio: " << *it << std::endl << "bailing..." << std::endl;
598  exit(1);
599  }
600  } else if (deck.hasKeyword("LAMEMOD") && deck.hasKeyword("SHEARMOD")) {
601  std::vector<double> lame = deck.getKeyword("LAMEMOD").getRawDoubleData();
602  std::vector<double> shear = deck.getKeyword("SHEARMOD").getRawDoubleData();
603  Emod.resize(lame.size());
604  Poiss.resize(lame.size());
605  for (size_t i=0;i<lame.size();++i) {
606  Emod[i] = shear[i]*(3*lame[i]+2*shear[i])/(lame[i]+shear[i]);
607  Poiss[i] = 0.5*lame[i]/(lame[i]+shear[i]);
608  }
609  std::vector<double>::const_iterator it = std::min_element(Poiss.begin(), Poiss.end());
610  if (*it < 0) {
611  std::cerr << "Auxetic material specified for cell " << it-Poiss.begin() << std::endl
612  << "Lamè modulus: " << lame[it-Poiss.begin()] << " Shearmodulus: " << shear[it-Poiss.begin()] << std::endl
613  << "Emod: "<< Emod[it-Poiss.begin()] << " Poisson's ratio: " << *it << std::endl << "bailing..." << std::endl;
614  exit(1);
615  }
616  } else if (deck.hasKeyword("BULKMOD") && deck.hasKeyword("SHEARMOD")) {
617  std::vector<double> bulk = deck.getKeyword("BULKMOD").getRawDoubleData();
618  std::vector<double> shear = deck.getKeyword("SHEARMOD").getRawDoubleData();
619  Emod.resize(bulk.size());
620  Poiss.resize(bulk.size());
621  for (size_t i=0;i<bulk.size();++i) {
622  Emod[i] = 9*bulk[i]*shear[i]/(3*bulk[i]+shear[i]);
623  Poiss[i] = 0.5*(3*bulk[i]-2*shear[i])/(3*bulk[i]+shear[i]);
624  }
625  std::vector<double>::const_iterator it = std::min_element(Poiss.begin(), Poiss.end());
626  if (*it < 0) {
627  std::cerr << "Auxetic material specified for cell " << it-Poiss.begin() << std::endl
628  << "Bulkmodulus: " << bulk[it-Poiss.begin()] << " Shearmodulus: " << shear[it-Poiss.begin()] << std::endl
629  << "Emod: "<< Emod[it-Poiss.begin()] << " Poisson's ratio: " << *it << std::endl << "bailing..." << std::endl;
630  exit(1);
631  }
632  } else if (deck.hasKeyword("PERMX") && deck.hasKeyword("PORO")) {
633  std::cerr << "WARNING: Using PERMX and PORO for elastic material properties" << std::endl;
634  Emod = deck.getKeyword("PERMX").getRawDoubleData();
635  Poiss = deck.getKeyword("PORO").getRawDoubleData();
636  } else {
637  std::cerr << "No material data found in eclipse file, aborting" << std::endl;
638  exit(1);
639  }
640  if (deck.hasKeyword("SATNUM"))
641  satnum = deck.getKeyword("SATNUM").getIntData();
642  if (deck.hasKeyword("RHO"))
643  rho = deck.getKeyword("RHO").getRawDoubleData();
644  }
645  // scale E modulus of materials
646  if (Escale > 0) {
647  Emin = *std::min_element(Emod.begin(),Emod.end());
648  for (size_t i=0;i<Emod.size();++i)
649  Emod[i] *= Escale/Emin;
650  }
651 
652  // make sure we only instance a minimal amount of materials.
653  // also map the correct material to the correct cells.
654  // their original ordering is as in the eclipse file itself
655  // while globalCell holds the map of cells kept after preprocessing
656  // the grid
657  std::vector<int> cells = gv.globalCell();
658  int j=0;
659  std::map<Material*,double> volume;
660  for (size_t i=0;i<cells.size();++i) {
661  int k = cells[i];
662  MaterialMap::iterator it;
663  if ((it = cache.find(std::make_pair(Emod[k],Poiss[k]))) != cache.end()) {
664  assert(gv.cellVolume(i) > 0);
665  volume[it->second.get()] += gv.cellVolume(i);
666  materials.push_back(it->second);
667  } else {
668  std::shared_ptr<Material> mat(new Isotropic(j++,Emod[k],Poiss[k]));
669  cache.insert(std::make_pair(std::make_pair(Emod[k],Poiss[k]),mat));
670  assert(gv.cellVolume(i) > 0);
671  volume.insert(std::make_pair(mat.get(),gv.cellVolume(i)));
672  materials.push_back(mat);
673  }
674  if (!satnum.empty()) {
675  if (satnum[k] > (int)volumeFractions.size())
676  volumeFractions.resize(satnum[k]);
677  volumeFractions[satnum[k]-1] += gv.cellVolume(i);
678  }
679  if (!rho.empty())
680  upscaledRho += gv.cellVolume(i)*rho[k];
681  }
682  std::cout << "Number of materials: " << cache.size() << std::endl;
683 
684  double totalvolume=0;
685  for (std::map<Material*,double>::iterator it = volume.begin();
686  it != volume.end(); ++it)
687  totalvolume += it->second;
688 
689  // statistics
690  if (satnum.empty()) {
691  int i=0;
692  for (MaterialMap::iterator it = cache.begin(); it != cache.end(); ++it, ++i) {
693  std::cout << " Material" << i+1 << ": " << 100.f*volume[it->second.get()]/totalvolume << '%' << std::endl;
694  volumeFractions.push_back(volume[it->second.get()]/totalvolume);
695  }
696  bySat = false;
697  } else {
698  for (size_t jj=0; jj < volumeFractions.size(); ++jj) {
699  volumeFractions[jj] /= totalvolume;
700  std::cout << "SATNUM " << jj+1 << ": " << volumeFractions[jj]*100 << '%' << std::endl;
701  }
702  bySat = true;
703  }
704  if (upscaledRho > 0) {
705  upscaledRho /= totalvolume;
706  std::cout << "Upscaled density: " << upscaledRho << std::endl;
707  }
708 }
709 
710 IMPL_FUNC(void, loadMaterialsFromRocklist(const std::string& file,
711  const std::string& rocklist))
712 {
713  std::vector< std::shared_ptr<Material> > cache;
714  // parse the rocklist
715  std::ifstream f;
716  f.open(rocklist.c_str());
717  int mats;
718  f >> mats;
719  for (int i=0;i<mats;++i) {
720  std::string matfile;
721  f >> matfile;
722  cache.push_back(std::shared_ptr<Material>(Material::create(i+1,matfile)));
723  }
724 
725  // scale E modulus of materials
726  if (Escale > 0) {
727  Emin=1e10;
728  for (size_t i=0;i<cache.size();++i)
729  Emin = std::min(Emin,((Isotropic*)cache[i].get())->getE());
730  for (size_t i=0;i<cache.size();++i) {
731  double lE = ((Isotropic*)cache[i].get())->getE();
732  ((Isotropic*)cache[i].get())->setE(lE*Escale/Emin);
733  }
734  }
735  std::vector<double> volume;
736  volume.resize(cache.size());
737  if (file == "uniform") {
738  if (cache.size() == 1) {
739  for (int i=0;i<gv.size(0);++i)
740  materials.push_back(cache[0]);
741  volume[0] = 1;
742  } else {
743  for (int i=0;i<gv.size(0);++i) {
744  materials.push_back(cache[i % cache.size()]);
745  volume[i % cache.size()] += gv.cellVolume(i);
746  }
747  }
748  } else {
749  Opm::Parser parser;
750  auto deck = parser.parseFile(file);
751  std::vector<int> satnum = deck.getKeyword("SATNUM").getIntData();
752  std::vector<int> cells = gv.globalCell();
753  for (size_t i=0;i<cells.size();++i) {
754  int k = cells[i];
755  if (satnum[k]-1 >= (int)cache.size()) {
756  std::cerr << "Material " << satnum[k] << " referenced but not available. Check your rocklist." << std::endl;
757  exit(1);
758  }
759  materials.push_back(cache[satnum[k]-1]);
760  volume[satnum[k]-1] += gv.cellVolume(i);
761  }
762  }
763  std::cout << "Number of materials: " << cache.size() << std::endl;
764  // statistics
765  double totalvolume = std::accumulate(volume.begin(),volume.end(),0.f);
766  for (size_t i=0;i<cache.size();++i) {
767  std::cout << " SATNUM " << i+1 << ": " << 100.f*volume[i]/totalvolume << '%' << std::endl;
768  volumeFractions.push_back(volume[i]/totalvolume);
769  }
770  bySat = true;
771 }
772 
773 IMPL_FUNC(void, fixCorners(const double* min, const double* max))
774 {
775  ctype c[8][3] = {{min[0],min[1],min[2]},
776  {max[0],min[1],min[2]},
777  {min[0],max[1],min[2]},
778  {max[0],max[1],min[2]},
779  {min[0],min[1],max[2]},
780  {max[0],min[1],max[2]},
781  {min[0],max[1],max[2]},
782  {max[0],max[1],max[2]}};
783  for (int i=0;i<8;++i) {
784  GlobalCoordinate coord;
785  coord[0] = c[i][0]; coord[1] = c[i][1]; coord[2] = c[i][2];
786  fixPoint(XYZ,coord);
787  }
788 }
789 
790 IMPL_FUNC(void, periodicBCs(const double* min, const double* max))
791 {
792  // this method
793  // 1. fixes the primal corner dofs
794  // 2. extracts establishes a quad grid for the left and front sides,
795  // while a point grid is created for the right and back sides.
796  // 3. establishes strong couplings (MPC)
797 
798  // step 1
799  fixCorners(min,max);
800 
801  // step 2
802  determineSideFaces(min,max);
803  std::cout << "Xslave " << slave[0].size() << " "
804  << "Yslave " << slave[1].size() << " "
805  << "Zslave " << slave[2].size() << std::endl;
806  std::cout << "Xmaster " << master[0].size() << " "
807  << "Ymaster " << master[1].size() << " "
808  << "Zmaster " << master[2].size() << std::endl;
809 
810  // step 3
811  periodicPlane(X,XYZ,slave[0],master[0]);
812  periodicPlane(Y,XYZ,slave[1],master[1]);
813  periodicPlane(Z,XYZ,slave[2],master[2]);
814 }
815 
816 IMPL_FUNC(void, periodicBCsMortar(const double* min,
817  const double* max,
818  int n1, int n2,
819  int p1, int p2))
820 {
821  // this method
822  // 1. fixes the primal corner dofs
823  // 2. establishes strong couplings (MPC) on top and bottom
824  // 3. extracts and establishes a quad grid for the left/right/front/back sides
825  // 4. establishes grids for the dual dofs
826  // 5. calculates the coupling matrix between the left/right sides
827  // 6. calculates the coupling matrix between the front/back sides
828  //
829  // The Mortar linear system is of the form
830  // [A B] [u] [b]
831  // [B' 0] [l] = [0]
832 
833  // step 1
834  fixCorners(min,max);
835 
836  // step 2
837  std::cout << "\textracting nodes on top face..." << std::endl;
838  slave.push_back(extractFace(Z,max[2]));
839  std::cout << "\treconstructing bottom face..." << std::endl;
840  BoundaryGrid bottom = extractMasterFace(Z,min[2]);
841  std::cout << "\testablishing couplings on top/bottom..." << std::endl;
842  periodicPlane(Z,XYZ,slave[0],bottom);
843  std::cout << "\tinitializing matrix..." << std::endl;
844  A.initForAssembly();
845 
846  // step 3
847  std::cout << "\treconstructing left face..." << std::endl;
848  master.push_back(extractMasterFace(X, min[0], LEFT, true));
849  std::cout << "\treconstructing right face..." << std::endl;
850  master.push_back(extractMasterFace(X, max[0], RIGHT, true));
851  std::cout << "\treconstructing front face..." << std::endl;
852  master.push_back(extractMasterFace(Y, min[1], LEFT, true));
853  std::cout << "\treconstructing back face..." << std::endl;
854  master.push_back(extractMasterFace(Y, max[1], RIGHT, true));
855 
856  std::cout << "\testablished YZ multiplier grid with " << n2 << "x1" << " elements" << std::endl;
857 
858  // step 4
859  BoundaryGrid::FaceCoord fmin,fmax;
860  fmin[0] = min[1]; fmin[1] = min[2];
861  fmax[0] = max[1]; fmax[1] = max[2];
862  BoundaryGrid lambdax = BoundaryGrid::uniform(fmin, fmax, n2, 1, true);
863 
864  fmin[0] = min[0]; fmin[1] = min[2];
865  fmax[0] = max[0]; fmax[1] = max[2];
866  std::cout << "\testablished XZ multiplier grid with " << n1 << "x1" << " elements" << std::endl;
867  BoundaryGrid lambday = BoundaryGrid::uniform(fmin, fmax, n1, 1, true);
868 
869  addBBlockMortar(master[0], lambdax, 0, 1, p2, 0);
870  int eqns = addBBlockMortar(master[1], lambdax, 0, 1, p2, 0);
871  addBBlockMortar(master[2], lambday, 1, 1, p1, eqns);
872  int eqns2 = addBBlockMortar(master[3], lambday, 1, 1, p1, eqns);
873 
874  MatrixOps::fromAdjacency(B, Bpatt, A.getEqns(), eqns+eqns2);
875  Bpatt.clear();
876 
877  // step 5
878  std::cout << "\tassembling YZ mortar matrix..." << std::endl;
879  assembleBBlockMortar(master[0], lambdax, 0, 1, p2, 0);
880  assembleBBlockMortar(master[1], lambdax, 0, 1, p2, 0, -1.0);
881 
882  // step 6
883  std::cout << "\tassembling XZ mortar matrix..." << std::endl;
884  assembleBBlockMortar(master[2], lambday, 1, 1, p1, eqns);
885  assembleBBlockMortar(master[3], lambday, 1, 1, p1, eqns, -1.0);
886 
887  master.clear();
888  slave.clear();
889 }
890 
891  template<class M, class A>
892 static void applyMortarBlock(int i, const Matrix& B, M& T,
893  A& upre)
894 {
895  Vector v, v2;
896  v.resize(B.N());
897  v2.resize(B.N());
898  v = 0;
899  v2 = 0;
900  MortarBlockEvaluator<Dune::Preconditioner<Vector,Vector> > pre(upre, B);
901  v[i] = 1;
902  pre.apply(v, v2);
903  for (size_t j=0; j < B.M(); ++j)
904  T[j][i] = v2[j];
905 }
906 
907 IMPL_FUNC(void, setupSolvers(const LinSolParams& params))
908 {
909  int siz = A.getOperator().N(); // system size
910  int numsolvers = 1;
911 #ifdef HAVE_OPENMP
912  numsolvers = omp_get_max_threads();
913 #endif
914 
915  if (params.type == ITERATIVE) {
916  op.reset(new Operator(A.getOperator()));
917  bool copy;
918  upre.push_back(PC::setup(params.steps[0], params.steps[1],
919  params.coarsen_target, params.zcells,
920  op, gv, A, copy));
921 
922  // Mortar in use
923  if (B.N()) {
924  siz += B.M();
925 
926  // schur system: B'*P^-1*B
927  if (params.mortarpre >= SCHUR) {
928  Dune::DynamicMatrix<double> T(B.M(), B.M());
929  std::cout << "\tBuilding preconditioner for multipliers..." << std::endl;
930  LoggerHelper help(B.M(), 5, 100);
931 
932  std::vector< std::shared_ptr<Dune::Preconditioner<Vector,Vector> > > pc;
933  pc.resize(B.M());
934 
935  if (params.mortarpre == SCHUR ||
936  (params.mortarpre == SCHURAMG &&
937  params.pre == AMG) ||
938  (params.mortarpre == SCHURSCHWARZ &&
939  params.pre == SCHWARZ) ||
940  (params.mortarpre == SCHURTWOLEVEL &&
941  params.pre == TWOLEVEL)) {
942  pc[0] = upre[0];
943  if (copy) {
944  for (size_t i=1;i<B.M();++i)
945  pc[i].reset(new PCType(*upre[0]));
946  }
947  } else if (params.mortarpre == SCHURAMG) {
948  std::shared_ptr<AMG1<SSORSmoother>::type> mpc;
949  pc[0] = mpc = AMG1<SSORSmoother>::setup(params.steps[0],
950  params.steps[1],
951  params.coarsen_target,
952  params.zcells,
953  op, gv, A, copy);
954  for (size_t i=1;i<B.M();++i)
955  pc[i].reset(new AMG1<SSORSmoother>::type(*mpc));
956  } else if (params.mortarpre == SCHURSCHWARZ) {
957  std::shared_ptr<Schwarz::type> mpc;
958  pc[0] = mpc = Schwarz::setup(params.steps[0],
959  params.steps[1],
960  params.coarsen_target,
961  params.zcells,
962  op, gv, A, copy);
963  } else if (params.mortarpre == SCHURTWOLEVEL) {
964  std::shared_ptr<AMG2Level<SSORSmoother>::type> mpc;
965  pc[0] = mpc = AMG2Level<SSORSmoother>::setup(params.steps[0],
966  params.steps[1],
967  params.coarsen_target,
968  params.zcells,
969  op, gv, A, copy);
970  for (size_t i=1;i<B.M();++i)
971  pc[i].reset(new AMG2Level<SSORSmoother>::type(*mpc));
972  }
973  for (int t=0;t<5;++t) {
974 #pragma omp parallel for schedule(static)
975  for (int i=help.group(t).first; i < help.group(t).second; ++i)
976  applyMortarBlock(i, B, T, *pc[copy?i:0]);
977 
978  help.log(help.group(t).second,
979  "\t\t... still processing ... multiplier ");
980  }
981  P = MatrixOps::fromDense(T);
982  } else if (params.mortarpre == SIMPLE) {
983  Matrix D = MatrixOps::diagonal(A.getEqns());
984 
985  // scale by row sums
986  size_t row=0;
987  for (Matrix::ConstRowIterator it = A.getOperator().begin();
988  it != A.getOperator().end(); ++it, ++row) {
989  double alpha=0;
990  for (Matrix::ConstColIterator it2 = it->begin();
991  it2 != it->end(); ++it2) {
992  if (it2.index() != row)
993  alpha += fabs(*it2);
994  }
995  D[row][row] = 1.0/(A.getOperator()[row][row]/alpha);
996  }
997 
998  Matrix t1;
999  // t1 = Ad*B
1000  Dune::matMultMat(t1, D, B);
1001  // P = B'*t1 = B'*Ad*B
1002  Dune::transposeMatMultMat(P, B, t1);
1003  }
1004 
1005  if (params.uzawa) {
1006  std::shared_ptr<Dune::InverseOperator<Vector,Vector> > innersolver;
1007 
1008  innersolver.reset(new Dune::CGSolver<Vector>(*op, *upre[0], params.tol,
1009  params.maxit,
1010  verbose?2:(params.report?1:0)));
1011  op2.reset(new SchurEvaluator(*innersolver, B));
1012  lprep.reset(new LUSolver(P));
1013  lpre.reset(new SeqLU(*lprep));
1014  std::shared_ptr<Dune::InverseOperator<Vector,Vector> > outersolver;
1015  outersolver.reset(new Dune::CGSolver<Vector>(*op2, *lpre, params.tol*10,
1016  params.maxit,
1017  verbose?2:(params.report?1:0)));
1018  tsolver.push_back(SolverPtr(new UzawaSolver<Vector, Vector>(innersolver, outersolver, B)));
1019  } else {
1020  for (int i=0;i<numsolvers;++i) {
1021  if (copy && i != 0)
1022  upre.push_back(PCPtr(new PCType(*upre[0])));
1023  tmpre.push_back(MortarAmgPtr(new MortarSchurPre<PCType>(P, B,
1024  *upre[copy?i:0],
1025  params.symmetric)));
1026  }
1027  meval.reset(new MortarEvaluator(A.getOperator(), B));
1028  if (params.symmetric) {
1029  for (int i=0;i<numsolvers;++i)
1030  tsolver.push_back(SolverPtr(new Dune::MINRESSolver<Vector>(*meval, *tmpre[i],
1031  params.tol,
1032  params.maxit,
1033  verbose?2:(params.report?1:0))));
1034  } else {
1035  for (int i=0;i<numsolvers;++i)
1036  tsolver.push_back(SolverPtr(new Dune::RestartedGMResSolver<Vector>(*meval, *tmpre[i],
1037  params.tol,
1038  params.restart,
1039  params.maxit,
1040  verbose?2:(params.report?1:0))));
1041  }
1042  }
1043  } else {
1044  for (int i=0;i<numsolvers;++i) {
1045  if (copy && i != 0)
1046  upre.push_back(PCPtr(new PCType(*upre[0])));
1047  tsolver.push_back(SolverPtr(new Dune::CGSolver<Vector>(*op, *upre[copy?i:0],
1048  params.tol,
1049  params.maxit,
1050  verbose?2:(params.report?1:0))));
1051  }
1052  }
1053  } else {
1054  if (B.N())
1055  A.getOperator() = MatrixOps::augment(A.getOperator(), B,
1056  0, A.getOperator().M(), true);
1057 #if HAVE_UMFPACK || HAVE_SUPERLU
1058  tsolver.push_back(SolverPtr(new LUSolver(A.getOperator(),
1059  verbose?2:(params.report?1:0))));
1060  for (int i=1;i<numsolvers;++i)
1061  tsolver.push_back(tsolver.front());
1062 #else
1063  std::cerr << "No direct solver available" << std::endl;
1064  exit(1);
1065 #endif
1066  siz = A.getOperator().N();
1067  }
1068 
1069  for (int i=0;i<6;++i)
1070  b[i].resize(siz);
1071 }
1072 
1073 IMPL_FUNC(void, solve(int loadcase))
1074 {
1075  try {
1076  Dune::InverseOperatorResult r;
1077 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
1078  u[loadcase].resize(b[loadcase].size());
1079 #else
1080  u[loadcase].resize(b[loadcase].size(), false);
1081 #endif
1082  u[loadcase] = 0;
1083  int solver=0;
1084 #ifdef HAVE_OPENMP
1085  solver = omp_get_thread_num();
1086 #endif
1087 
1088  tsolver[solver]->apply(u[loadcase], b[loadcase], r);
1089 
1090  std::cout << "\tsolution norm: " << u[loadcase].two_norm() << std::endl;
1091  } catch (Dune::ISTLError& e) {
1092  std::cerr << "exception thrown " << e << std::endl;
1093  }
1094 }
1095 
1096 }} // namespace Opm, Elasticity
1097 
1098 #endif
BoundaryGrid::Vertex maxXmaxY(std::vector< BoundaryGrid::Vertex > &in)
Find the vertex in the vector with maximum X and maximum Y.
Definition: boundarygrid.cpp:371
BoundaryGrid::Vertex minXminY(std::vector< BoundaryGrid::Vertex > &in)
Find the vertex in the vector with minimum X and minimum Y.
Definition: boundarygrid.cpp:351
BoundaryGrid::Vertex maxXminY(std::vector< BoundaryGrid::Vertex > &in)
Find the vertex in the vector with maximum X and minimum Y.
Definition: boundarygrid.cpp:361
BoundaryGrid::Vertex minXmaxY(std::vector< BoundaryGrid::Vertex > &in)
Find the vertex in the vector with minimum X and maximum Y.
Definition: boundarygrid.cpp:381
Helper class for progress logging during time consuming processes.
Definition: logutils.hpp:21
static void extract(FaceCoord &res, const GlobalCoordinate &coord, int dir)
Helper function for extracting given 2D coordinates from a 3D coordinate.
Definition: boundarygrid.cpp:70
Dune::FieldVector< double, 2 > FaceCoord
A coordinate on the quad grid.
Definition: boundarygrid.hh:38
static BoundaryGrid uniform(const FaceCoord &min, const FaceCoord &max, int k1, int k2, bool dc=false)
Establish an uniform quad grid.
Definition: boundarygrid.cpp:24
static Material * create(int ID, const Dune::DynamicVector< double > &params)
Creates a material object of a given type.
Definition: material.cpp:23
static Matrix fromDense(const Dune::DynamicMatrix< double > &T)
Create a sparse matrix from a dense matrix.
Definition: matrixops.cpp:49
static void fromAdjacency(Matrix &A, const AdjacencyPattern &adj, int rows, int cols)
Create a sparse matrix from a given adjacency pattern.
Definition: matrixops.cpp:25
static Matrix augment(const Matrix &A, const Matrix &B, size_t r0, size_t c0, bool symmetric)
Augment a matrix with another.
Definition: matrixops.cpp:126
static Matrix diagonal(size_t N)
Returns a diagonal matrix.
Definition: matrixops.cpp:198
static const P1ShapeFunctionSet & instance()
Get the only instance of this class.
Definition: shapefunctions.hpp:193
Dune::MatrixAdapter< Matrix, Vector, Vector > Operator
A linear operator.
Definition: elasticity_preconditioners.hpp:45
Dune::BCRSMatrix< Dune::FieldMatrix< double, 1, 1 > > Matrix
A sparse matrix holding our operator.
Definition: matrixops.hpp:27
Dune::BlockVector< Dune::FieldVector< double, 1 > > Vector
A vector holding our RHS.
Definition: matrixops.hpp:33
Inverting small matrices.
Definition: ImplicitAssembly.hpp:43
static std::shared_ptr< type > setup(int pre, int post, int target, int zcells, std::shared_ptr< Operator > &op, const Dune::CpGrid &, ASMHandler< Dune::CpGrid > &, bool &copy)
Setup preconditioner.
Definition: elasticity_preconditioners.hpp:121
static std::shared_ptr< type > setup(int pre, int post, int target, int zcells, std::shared_ptr< Operator > &op, const Dune::CpGrid &gv, ASMHandler< Dune::CpGrid > &A, bool &copy)
Setup preconditioner.
Definition: elasticity_preconditioners.hpp:184
static std::shared_ptr< type > setup(int, int, int, int, std::shared_ptr< Operator > &op, const Dune::CpGrid &gv, ASMHandler< Dune::CpGrid > &A, bool &copy)
Setup preconditioner.
Definition: elasticity_preconditioners.hpp:82