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