12 #ifndef OPM_ELASTICITY_UPSCALE_IMPL_HPP
13 #define OPM_ELASTICITY_UPSCALE_IMPL_HPP
21 #include <opm/parser/eclipse/Deck/DeckKeyword.hpp>
22 #include <opm/parser/eclipse/Parser/ParseContext.hpp>
25 namespace Elasticity {
28 #define IMPL_FUNC(A,B) template<class GridType, class PC> \
29 A ElasticityUpscale<GridType, PC>::B
31 IMPL_FUNC(std::vector<BoundaryGrid::Vertex>,
32 extractFace(Direction dir, ctype coord))
34 std::vector<BoundaryGrid::Vertex> result;
35 const LeafVertexIterator itend = gv.leafGridView().template end<dim>();
38 Dune::LeafMultipleCodimMultipleGeomTypeMapper<GridType,
39 Dune::MCMGVertexLayout> mapper(gv);
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);
55 IMPL_FUNC(BoundaryGrid, extractMasterFace(Direction dir,
60 static const int V1[3][4] = {{0,2,4,6},
63 static const int V2[3][4] = {{1,3,5,7},
66 const LeafIndexSet& set = gv.leafGridView().indexSet();
69 int i = log2(
float(dir));
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;
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) {
86 idx = set.subIndex(*cell,V1[i][j],dim);
88 idx = set.subIndex(*cell,V2[i][j],dim);
89 pos = gv.vertexPosition(idx);
90 if (!isOnPlane(dir,pos,coord))
92 BoundaryGrid::Vertex v;
98 if (verts.size() == 4) {
100 q.v[0] = minXminY(verts);
101 q.v[1] = maxXminY(verts);
103 q.v[2] = minXmaxY(verts);
104 q.v[3] = maxXmaxY(verts);
106 q.v[2] = maxXmaxY(verts);
107 q.v[3] = minXmaxY(verts);
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);
116 if (it == nodeMap.end())
117 nodeMap[q.v[0].c[0]].push_back(q);
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]);
133 IMPL_FUNC(
void, determineSideFaces(
const double* min,
const double* max))
135 master.push_back(extractMasterFace(X,min[0]));
136 master.push_back(extractMasterFace(Y,min[1]));
137 master.push_back(extractMasterFace(Z,min[2]));
139 slave.push_back(extractFace(X,max[0]));
140 slave.push_back(extractFace(Y,max[1]));
141 slave.push_back(extractFace(Z,max[2]));
144 IMPL_FUNC(
void, findBoundaries(
double* min,
double* max))
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>();
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]);
160 IMPL_FUNC(
void, addMPC(Direction dir,
int slavenode,
161 const BoundaryGrid::Vertex& m))
163 MPC* mpc =
new MPC(slavenode,log2(
float(dir))+1);
165 mpc->addMaster(m.i,log2(
float(dir))+1,1.f);
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]);
174 IMPL_FUNC(
void, periodicPlane(Direction ,
176 const std::vector<BoundaryGrid::Vertex>& slavepointgrid,
177 const BoundaryGrid& mastergrid))
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);
194 static std::vector< std::vector<int> > renumber(
const BoundaryGrid& b,
198 std::vector<std::vector<int> > nodes;
199 nodes.resize(b.size());
205 nodes[0].push_back(-1);
206 for (
size_t e=0; e < b.size(); ++e) {
208 for (
int i2=0; i2 < n2; ++i2) {
210 nodes[e].push_back(nodes[e-1][i2*n1+n1-1]);
212 int start = (e==0 && i2 != 0)?0:1;
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]);
220 for (
int i1=start; i1 < n1; ++i1) {
222 nodes[e].push_back(nodes[0][i2*n1]);
224 nodes[e].push_back(totalDOFs++);
233 IMPL_FUNC(
int, addBBlockMortar(
const BoundaryGrid& b1,
234 const BoundaryGrid& interface,
235 int ,
int n1,
int n2,
240 std::vector<std::vector<int> > lnodes = renumber(interface, n1+1,
243 Bpatt.resize(A.getEqns());
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);
252 for (
size_t n=0;n<mpc->getNoMaster();++n) {
253 int dof = A.getEquationForDof(mpc->getMaster(n).node,d);
255 for (
size_t j=0;j<lnodes[p].size();++j) {
256 int indexj = 3*lnodes[p][j]+d;
258 Bpatt[dof].insert(indexj+colofs);
263 int dof = A.getEquationForDof(b1.getQuad(p,q).v[i].i,d);
265 for (
size_t j=0;j<lnodes[p].size();++j) {
266 int indexj = 3*lnodes[p][j]+d;
268 Bpatt[dof].insert(indexj+colofs);
280 IMPL_FUNC(
void, assembleBBlockMortar(
const BoundaryGrid& b1,
281 const BoundaryGrid& interface,
287 P1ShapeFunctionSet<ctype,ctype,2> ubasis =
291 PNShapeFunctionSet<2> lbasis(n1+1, n2+1);
295 std::vector<std::vector<int> > lnodes = renumber(interface, n1+1,
298 Dune::GeometryType gt;
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);
307 typename Dune::QuadratureRule<ctype,2>::const_iterator r;
308 Dune::DynamicMatrix<ctype> E(ubasis.n,(n1+1)*(n2+1),0.0);
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);
318 for (r = rule.begin(); r != rule.end();++r) {
319 ctype detJ = hex.integrationElement(r->position());
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();
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);
339 for (
size_t n=0;n<mpc->getNoMaster();++n) {
340 int indexi = A.getEquationForDof(mpc->getMaster(n).node,d);
342 for (
size_t j=0;j<lnodes[p].size();++j) {
343 int indexj = lnodes[p][j]*3+d;
345 B[indexi][indexj+colofs] += alpha*E[i][j];
350 int indexi = A.getEquationForDof(qu.v[i].i,d);
352 for (
size_t j=0;j<lnodes[p].size();++j) {
353 int indexj = lnodes[p][j]*3+d;
355 B[indexi][indexj+colofs] += alpha*E[i][j];
363 help.log(g,
"\t\t\t... still processing ... pillar ");
367 IMPL_FUNC(
void, fixPoint(Direction dir,
368 GlobalCoordinate coord,
369 const NodeValue& value))
371 typedef typename GridType::LeafGridView::template Codim<dim>::Iterator VertexLeafIterator;
372 const VertexLeafIterator itend = gv.leafGridView().template end<dim>();
375 Dune::LeafMultipleCodimMultipleGeomTypeMapper<GridType,
376 Dune::MCMGVertexLayout> mapper(gv);
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));
387 IMPL_FUNC(
bool, isOnPlane(Direction plane,
388 GlobalCoordinate coord,
391 if (plane < X || plane > Z)
393 int p = log2(
float(plane));
394 ctype delta = fabs(value-coord[p]);
398 IMPL_FUNC(
void, fixLine(Direction dir,
400 const NodeValue& value))
402 typedef typename GridType::LeafGridView::template Codim<dim>::Iterator VertexLeafIterator;
403 const VertexLeafIterator itend = gv.leafGridView().template end<dim>();
406 Dune::LeafMultipleCodimMultipleGeomTypeMapper<GridType,
407 Dune::MCMGVertexLayout> mapper(gv);
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));
418 IMPL_FUNC(
bool, isOnLine(Direction dir,
419 GlobalCoordinate coord,
422 if (dir < X || dir > Z)
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)
430 if (delta > tol || delta < -tol)
436 IMPL_FUNC(
bool, isOnPoint(GlobalCoordinate coord,
437 GlobalCoordinate point))
439 GlobalCoordinate delta = point-coord;
440 return delta.one_norm() < tol;
443 IMPL_FUNC(
void, assemble(
int loadcase,
bool matrix))
445 const int comp = 3+(dim-2)*3;
446 static const int bfunc = 4+(dim-2)*4;
448 Dune::FieldVector<ctype,comp> eps0;
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;
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)
476 materials[color[i][j][k]]->getConstitutiveMatrix(C);
478 Dune::GeometryType gt = it->type();
480 Dune::FieldMatrix<ctype,dim*bfunc,dim*bfunc> Aq;
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) {
489 Dune::FieldMatrix<ctype,dim,dim> jacInvTra =
490 it->geometry().jacobianInverseTransposed(r->position());
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;
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;
501 Dune::FieldMatrix<ctype,comp,dim*bfunc> B;
502 E.getBmatrix(B,r->position(),jacInvTra);
505 E.getStiffnessMatrix(Aq,B,C,detJ*r->weight());
511 Dune::FieldVector<ctype,dim*bfunc> temp;
512 temp = Dune::FMatrixHelp::multTransposed(B,Dune::FMatrixHelp::mult(C,eps0));
513 temp *= -detJ*r->weight();
517 A.addElement(KP,EP,it,(loadcase > -1)?&b[loadcase]:NULL);
523 IMPL_FUNC(template<int comp>
void,
524 averageStress(Dune::FieldVector<ctype,comp>& sigma,
525 const Vector& uarg,
int loadcase))
527 if (loadcase < 0 || loadcase > 5)
530 static const int bfunc = 4+(dim-2)*4;
532 const LeafIterator itend = gv.leafGridView().template end<0>();
534 Dune::FieldMatrix<ctype,comp,comp> C;
535 Dune::FieldVector<ctype,comp> eps0;
541 for (LeafIterator it = gv.leafGridView().template begin<0>(); it != itend; ++it) {
542 materials[m++]->getConstitutiveMatrix(C);
544 Dune::GeometryType gt = it->type();
546 Dune::FieldVector<ctype,bfunc*dim> v;
547 A.extractValues(v,uarg,it);
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) {
554 Dune::FieldMatrix<ctype,dim,dim> jacInvTra =
555 it->geometry().jacobianInverseTransposed(r->position());
557 ctype detJ = it->geometry().integrationElement(r->position());
559 volume += detJ*r->weight();
561 Dune::FieldMatrix<ctype,comp,dim*bfunc> B;
562 E.getBmatrix(B,r->position(),jacInvTra);
564 Dune::FieldVector<ctype,comp> s;
565 E.getStressVector(s,v,eps0,B,C);
566 s *= detJ*r->weight();
572 sigma /= Escale/Emin;
575 IMPL_FUNC(
void, loadMaterialsFromGrid(
const std::string& file))
577 typedef std::map<std::pair<double,double>,
578 std::shared_ptr<Material> > MaterialMap;
580 std::vector<double> Emod;
581 std::vector<double> Poiss;
582 std::vector<int> satnum;
583 std::vector<double> rho;
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);
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());
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;
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]);
612 std::vector<double>::const_iterator it = std::min_element(Poiss.begin(), Poiss.end());
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;
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]);
628 std::vector<double>::const_iterator it = std::min_element(Poiss.begin(), Poiss.end());
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;
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();
640 std::cerr <<
"No material data found in eclipse file, aborting" << std::endl;
643 if (deck.hasKeyword(
"SATNUM"))
644 satnum = deck.getKeyword(
"SATNUM").getIntData();
645 if (deck.hasKeyword(
"RHO"))
646 rho = deck.getKeyword(
"RHO").getRawDoubleData();
650 Emin = *std::min_element(Emod.begin(),Emod.end());
651 for (
size_t i=0;i<Emod.size();++i)
652 Emod[i] *= Escale/Emin;
660 std::vector<int> cells = gv.globalCell();
662 std::map<Material*,double> volume;
663 for (
size_t i=0;i<cells.size();++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);
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);
677 if (!satnum.empty()) {
678 if (satnum[k] > (
int)volumeFractions.size())
679 volumeFractions.resize(satnum[k]);
680 volumeFractions[satnum[k]-1] += gv.cellVolume(i);
683 upscaledRho += gv.cellVolume(i)*rho[k];
685 std::cout <<
"Number of materials: " << cache.size() << std::endl;
687 double totalvolume=0;
688 for (std::map<Material*,double>::iterator it = volume.begin();
689 it != volume.end(); ++it)
690 totalvolume += it->second;
693 if (satnum.empty()) {
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);
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;
707 if (upscaledRho > 0) {
708 upscaledRho /= totalvolume;
709 std::cout <<
"Upscaled density: " << upscaledRho << std::endl;
713 IMPL_FUNC(
void, loadMaterialsFromRocklist(
const std::string& file,
714 const std::string& rocklist))
716 std::vector< std::shared_ptr<Material> > cache;
719 f.open(rocklist.c_str());
722 for (
int i=0;i<mats;++i) {
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);
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]);
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);
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) {
759 if (satnum[k]-1 >= (
int)cache.size()) {
760 std::cerr <<
"Material " << satnum[k] <<
" referenced but not available. Check your rocklist." << std::endl;
763 materials.push_back(cache[satnum[k]-1]);
764 volume[satnum[k]-1] += gv.cellVolume(i);
767 std::cout <<
"Number of materials: " << cache.size() << std::endl;
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);
777 IMPL_FUNC(
void, fixCorners(
const double* min,
const double* max))
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];
794 IMPL_FUNC(
void, periodicBCs(
const double* min,
const double* max))
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;
815 periodicPlane(X,XYZ,slave[0],master[0]);
816 periodicPlane(Y,XYZ,slave[1],master[1]);
817 periodicPlane(Z,XYZ,slave[2],master[2]);
820 IMPL_FUNC(
void, periodicBCsMortar(
const double* min,
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;
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));
860 std::cout <<
"\testablished YZ multiplier grid with " << n2 <<
"x1" <<
" elements" << std::endl;
864 fmin[0] = min[1]; fmin[1] = min[2];
865 fmax[0] = max[1]; fmax[1] = max[2];
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;
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);
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);
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);
895 template<
class M,
class A>
896 static void applyMortarBlock(
int i,
const Matrix& B, M& T,
904 MortarBlockEvaluator<Dune::Preconditioner<Vector,Vector> > pre(upre, B);
907 for (
size_t j=0; j < B.M(); ++j)
911 IMPL_FUNC(
void, setupSolvers(
const LinSolParams& params))
913 int siz = A.getOperator().N();
916 numsolvers = omp_get_max_threads();
919 if (params.type == ITERATIVE) {
920 op.reset(
new Operator(A.getOperator()));
922 upre.push_back(PC::setup(params.steps[0], params.steps[1],
923 params.coarsen_target, params.zcells,
931 if (params.mortarpre >= SCHUR) {
932 Dune::DynamicMatrix<double> T(B.M(), B.M());
933 std::cout <<
"\tBuilding preconditioner for multipliers..." << std::endl;
936 std::vector< std::shared_ptr<Dune::Preconditioner<Vector,Vector> > > pc;
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)) {
948 for (
size_t i=1;i<B.M();++i)
949 pc[i].reset(
new PCType(*upre[0]));
951 }
else if (params.mortarpre == SCHURAMG) {
952 std::shared_ptr<AMG1<SSORSmoother>::type> mpc;
955 params.coarsen_target,
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;
964 params.coarsen_target,
967 }
else if (params.mortarpre == SCHURTWOLEVEL) {
968 std::shared_ptr<AMG2Level<SSORSmoother>::type> mpc;
971 params.coarsen_target,
974 for (
size_t i=1;i<B.M();++i)
975 pc[i].reset(
new AMG2Level<SSORSmoother>::type(*mpc));
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]);
982 help.log(help.group(t).second,
983 "\t\t... still processing ... multiplier ");
986 }
else if (params.mortarpre == SIMPLE) {
991 for (Matrix::ConstRowIterator it = A.getOperator().begin();
992 it != A.getOperator().end(); ++it, ++row) {
994 for (Matrix::ConstColIterator it2 = it->begin();
995 it2 != it->end(); ++it2) {
996 if (it2.index() != row)
999 D[row][row] = 1.0/(A.getOperator()[row][row]/alpha);
1004 Dune::matMultMat(t1, D, B);
1006 Dune::transposeMatMultMat(P, B, t1);
1010 std::shared_ptr<Dune::InverseOperator<Vector,Vector> > innersolver;
1012 innersolver.reset(
new Dune::CGSolver<Vector>(*op, *upre[0], params.tol,
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,
1021 verbose?2:(params.report?1:0)));
1022 tsolver.push_back(SolverPtr(
new UzawaSolver<Vector, Vector>(innersolver, outersolver, B)));
1024 for (
int i=0;i<numsolvers;++i) {
1026 upre.push_back(PCPtr(
new PCType(*upre[0])));
1027 tmpre.push_back(MortarAmgPtr(
new MortarSchurPre<PCType>(P, B,
1029 params.symmetric)));
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],
1037 verbose?2:(params.report?1:0))));
1039 for (
int i=0;i<numsolvers;++i)
1040 tsolver.push_back(SolverPtr(
new Dune::RestartedGMResSolver<Vector>(*meval, *tmpre[i],
1044 verbose?2:(params.report?1:0))));
1048 for (
int i=0;i<numsolvers;++i) {
1050 upre.push_back(PCPtr(
new PCType(*upre[0])));
1051 tsolver.push_back(SolverPtr(
new Dune::CGSolver<Vector>(*op, *upre[copy?i:0],
1054 verbose?2:(params.report?1:0))));
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());
1067 std::cerr <<
"No direct solver available" << std::endl;
1070 siz = A.getOperator().N();
1073 for (
int i=0;i<6;++i)
1077 IMPL_FUNC(
void, solve(
int loadcase))
1080 Dune::InverseOperatorResult r;
1081 u[loadcase].resize(b[loadcase].size(),
false);
1085 solver = omp_get_thread_num();
1088 tsolver[solver]->apply(u[loadcase], b[loadcase], r);
1090 std::cout <<
"\tsolution norm: " << u[loadcase].two_norm() << std::endl;
1091 }
catch (Dune::ISTLError& e) {
1092 std::cerr <<
"exception thrown " << e << std::endl;
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 ©)
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 ©)
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 > ¶ms)
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 ©)
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