12#ifndef OPM_ELASTICITY_UPSCALE_IMPL_HPP
13#define OPM_ELASTICITY_UPSCALE_IMPL_HPP
21#include <opm/input/eclipse/Deck/DeckKeyword.hpp>
27#define IMPL_FUNC(A,B) template<class GridType, class PC> \
28 A ElasticityUpscale<GridType, PC>::B
30IMPL_FUNC(std::vector<BoundaryGrid::Vertex>,
31 extractFace(Direction dir, ctype coord))
33 std::vector<BoundaryGrid::Vertex> result;
34 const LeafVertexIterator itend = gv.leafGridView().template end<dim>();
37 using LeafGridView = Dune::GridView<Dune::DefaultLeafGridViewTraits<GridType>>;
38 Dune::MultipleCodimMultipleGeomTypeMapper<LeafGridView> mapper(gv.leafGridView(), Dune::mcmgVertexLayout());
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);
54IMPL_FUNC(BoundaryGrid, extractMasterFace(Direction dir,
59 static const int V1[3][4] = {{0,2,4,6},
62 static const int V2[3][4] = {{1,3,5,7},
65 const LeafIndexSet& set = gv.leafGridView().indexSet();
68 int i = log2(
float(dir));
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;
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) {
85 idx = set.subIndex(*cell,V1[i][j],dim);
87 idx = set.subIndex(*cell,V2[i][j],dim);
88 pos = gv.vertexPosition(idx);
89 if (!isOnPlane(dir,pos,coord))
91 BoundaryGrid::Vertex v;
97 if (verts.size() == 4) {
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);
115 if (it == nodeMap.end())
116 nodeMap[q.v[0].c[0]].push_back(q);
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]);
132IMPL_FUNC(
void, determineSideFaces(
const double* min,
const double* max))
134 master.push_back(extractMasterFace(X,min[0]));
135 master.push_back(extractMasterFace(Y,min[1]));
136 master.push_back(extractMasterFace(Z,min[2]));
138 slave.push_back(extractFace(X,max[0]));
139 slave.push_back(extractFace(Y,max[1]));
140 slave.push_back(extractFace(Z,max[2]));
143IMPL_FUNC(
void, findBoundaries(
double* min,
double* max))
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>();
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]);
159IMPL_FUNC(
void, addMPC(Direction dir,
int slavenode,
160 const BoundaryGrid::Vertex& m))
162 MPC* mpc =
new MPC(slavenode,log2(
float(dir))+1);
164 mpc->addMaster(m.i,log2(
float(dir))+1,1.f);
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]);
173IMPL_FUNC(
void, periodicPlane(Direction ,
175 const std::vector<BoundaryGrid::Vertex>& slavepointgrid,
176 const BoundaryGrid& mastergrid))
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);
193static std::vector< std::vector<int> > renumber(
const BoundaryGrid& b,
197 std::vector<std::vector<int> > nodes;
198 nodes.resize(b.size());
204 nodes[0].push_back(-1);
205 for (
size_t e=0; e < b.size(); ++e) {
207 for (
int i2=0; i2 < n2; ++i2) {
209 nodes[e].push_back(nodes[e-1][i2*n1+n1-1]);
211 int start = (e==0 && i2 != 0)?0:1;
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]);
219 for (
int i1=start; i1 < n1; ++i1) {
221 nodes[e].push_back(nodes[0][i2*n1]);
223 nodes[e].push_back(totalDOFs++);
232IMPL_FUNC(
int, addBBlockMortar(
const BoundaryGrid& b1,
233 const BoundaryGrid& interface,
234 int ,
int n1,
int n2,
239 std::vector<std::vector<int> > lnodes = renumber(interface, n1+1,
242 Bpatt.resize(A.getEqns());
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);
251 for (
size_t n=0;n<mpc->getNoMaster();++n) {
252 int dof = A.getEquationForDof(mpc->getMaster(n).node,d);
254 for (
size_t j=0;j<lnodes[p].size();++j) {
255 int indexj = 3*lnodes[p][j]+d;
257 Bpatt[dof].insert(indexj+colofs);
262 int dof = A.getEquationForDof(b1.getQuad(p,q).v[i].i,d);
264 for (
size_t j=0;j<lnodes[p].size();++j) {
265 int indexj = 3*lnodes[p][j]+d;
267 Bpatt[dof].insert(indexj+colofs);
279IMPL_FUNC(
void, assembleBBlockMortar(
const BoundaryGrid& b1,
280 const BoundaryGrid& interface,
286 P1ShapeFunctionSet<ctype,ctype,2> ubasis =
290 PNShapeFunctionSet<2> lbasis(n1+1, n2+1);
294 std::vector<std::vector<int> > lnodes = renumber(interface, n1+1,
297 auto gt = Dune::GeometryTypes::cube(2);
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);
305 typename Dune::QuadratureRule<ctype,2>::const_iterator r;
306 Dune::DynamicMatrix<ctype> lE(ubasis.n,(n1+1)*(n2+1),0.0);
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);
316 for (r = rule.begin(); r != rule.end();++r) {
317 ctype detJ = hex.integrationElement(r->position());
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();
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);
337 for (
size_t n=0;n<mpc->getNoMaster();++n) {
338 int indexi = A.getEquationForDof(mpc->getMaster(n).node,d);
340 for (
size_t j=0;j<lnodes[p].size();++j) {
341 int indexj = lnodes[p][j]*3+d;
343 B[indexi][indexj+colofs] += alpha*lE[i][j];
348 int indexi = A.getEquationForDof(qu.v[i].i,d);
350 for (
size_t j=0;j<lnodes[p].size();++j) {
351 int indexj = lnodes[p][j]*3+d;
353 B[indexi][indexj+colofs] += alpha*lE[i][j];
361 help.log(g,
"\t\t\t... still processing ... pillar ");
365IMPL_FUNC(
void, fixPoint(Direction dir,
366 GlobalCoordinate coord,
367 const NodeValue& value))
369 typedef typename GridType::LeafGridView::template Codim<dim>::Iterator VertexLeafIterator;
370 const VertexLeafIterator itend = gv.leafGridView().template end<dim>();
373 using LeafGridView = Dune::GridView<Dune::DefaultLeafGridViewTraits<GridType>>;
374 Dune::MultipleCodimMultipleGeomTypeMapper<LeafGridView> mapper(gv.leafGridView(), Dune::mcmgVertexLayout());
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));
385IMPL_FUNC(
bool, isOnPlane(Direction plane,
386 GlobalCoordinate coord,
389 if (plane < X || plane > Z)
391 int p = log2(
float(plane));
392 ctype delta = fabs(value-coord[p]);
396IMPL_FUNC(
void, fixLine(Direction dir,
398 const NodeValue& value))
400 typedef typename GridType::LeafGridView::template Codim<dim>::Iterator VertexLeafIterator;
401 const VertexLeafIterator itend = gv.leafGridView().template end<dim>();
404 using LeafGridView = Dune::GridView<Dune::DefaultLeafGridViewTraits<GridType>>;
405 Dune::MultipleCodimMultipleGeomTypeMapper<LeafGridView> mapper(gv, Dune::mcmgVertexLayout());
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));
416IMPL_FUNC(
bool, isOnLine(Direction dir,
417 GlobalCoordinate coord,
420 if (dir < X || dir > Z)
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)
428 if (delta > tol || delta < -tol)
434IMPL_FUNC(
bool, isOnPoint(GlobalCoordinate coord,
435 GlobalCoordinate point))
437 GlobalCoordinate delta = point-coord;
438 return delta.one_norm() < tol;
441IMPL_FUNC(
void, assemble(
int loadcase,
bool matrix))
443 const int comp = 3+(dim-2)*3;
444 static const int bfunc = 4+(dim-2)*4;
446 Dune::FieldVector<ctype,comp> eps0;
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;
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)
474 materials[color[i][j][k]]->getConstitutiveMatrix(C);
476 Dune::GeometryType gt = it->type();
478 Dune::FieldMatrix<ctype,dim*bfunc,dim*bfunc> Aq;
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) {
487 Dune::FieldMatrix<ctype,dim,dim> jacInvTra =
488 it->geometry().jacobianInverseTransposed(r->position());
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;
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;
499 Dune::FieldMatrix<ctype,comp,dim*bfunc> lB;
500 E.getBmatrix(lB,r->position(),jacInvTra);
503 E.getStiffnessMatrix(Aq,lB,C,detJ*r->weight());
509 Dune::FieldVector<ctype,dim*bfunc> temp;
510 temp = Dune::FMatrixHelp::multTransposed(lB,Dune::FMatrixHelp::mult(C,eps0));
511 temp *= -detJ*r->weight();
515 A.addElement(KP,EP,it,(loadcase > -1)?&b[loadcase]:NULL);
521IMPL_FUNC(template<int comp>
void,
522 averageStress(Dune::FieldVector<ctype,comp>& sigma,
523 const Vector& uarg,
int loadcase))
525 if (loadcase < 0 || loadcase > 5)
528 static const int bfunc = 4+(dim-2)*4;
530 const LeafIterator itend = gv.leafGridView().template end<0>();
532 Dune::FieldMatrix<ctype,comp,comp> C;
533 Dune::FieldVector<ctype,comp> eps0;
539 for (LeafIterator it = gv.leafGridView().template begin<0>(); it != itend; ++it) {
540 materials[m++]->getConstitutiveMatrix(C);
542 Dune::GeometryType gt = it->type();
544 Dune::FieldVector<ctype,bfunc*dim> v;
545 A.extractValues(v,uarg,it);
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) {
552 Dune::FieldMatrix<ctype,dim,dim> jacInvTra =
553 it->geometry().jacobianInverseTransposed(r->position());
555 ctype detJ = it->geometry().integrationElement(r->position());
557 volume += detJ*r->weight();
559 Dune::FieldMatrix<ctype,comp,dim*bfunc> lB;
560 E.getBmatrix(lB,r->position(),jacInvTra);
562 Dune::FieldVector<ctype,comp> s;
563 E.getStressVector(s,v,eps0,lB,C);
564 s *= detJ*r->weight();
570 sigma /= Escale/Emin;
573IMPL_FUNC(
void, loadMaterialsFromGrid(
const std::string& file))
575 typedef std::map<std::pair<double,double>,
576 std::shared_ptr<Material> > MaterialMap;
578 std::vector<double> Emod;
579 std::vector<double> Poiss;
580 std::vector<int> satnum;
581 std::vector<double> rho;
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);
590 auto deck = parser.parseFile(file);
591 if (deck.hasKeyword(
"YOUNGMOD") && deck.hasKeyword(
"POISSONMOD")) {
592 Emod = deck[
"YOUNGMOD"].back().getRawDoubleData();
593 Poiss = deck[
"POISSONMOD"].back().getRawDoubleData();
594 std::vector<double>::const_iterator it = std::min_element(Poiss.begin(), Poiss.end());
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;
600 }
else if (deck.hasKeyword(
"LAMEMOD") && deck.hasKeyword(
"SHEARMOD")) {
601 std::vector<double> lame = deck[
"LAMEMOD"].back().getRawDoubleData();
602 std::vector<double> shear = deck[
"SHEARMOD"].back().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]);
609 std::vector<double>::const_iterator it = std::min_element(Poiss.begin(), Poiss.end());
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;
616 }
else if (deck.hasKeyword(
"BULKMOD") && deck.hasKeyword(
"SHEARMOD")) {
617 std::vector<double> bulk = deck[
"BULKMOD"].back().getRawDoubleData();
618 std::vector<double> shear = deck[
"SHEARMOD"].back().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]);
625 std::vector<double>::const_iterator it = std::min_element(Poiss.begin(), Poiss.end());
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;
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[
"PERMX"].back().getRawDoubleData();
635 Poiss = deck[
"PORO"].back().getRawDoubleData();
637 std::cerr <<
"No material data found in eclipse file, aborting" << std::endl;
640 if (deck.hasKeyword(
"SATNUM"))
641 satnum = deck[
"SATNUM"].back().getIntData();
642 if (deck.hasKeyword(
"RHO"))
643 rho = deck[
"RHO"].back().getRawDoubleData();
647 Emin = *std::min_element(Emod.begin(),Emod.end());
648 for (
size_t i=0;i<Emod.size();++i)
649 Emod[i] *= Escale/Emin;
657 std::vector<int> cells = gv.globalCell();
659 std::map<Material*,double> volume;
660 for (
size_t i=0;i<cells.size();++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);
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);
674 if (!satnum.empty()) {
675 if (satnum[k] > (
int)volumeFractions.size())
676 volumeFractions.resize(satnum[k]);
677 volumeFractions[satnum[k]-1] += gv.cellVolume(i);
680 upscaledRho += gv.cellVolume(i)*rho[k];
682 std::cout <<
"Number of materials: " << cache.size() << std::endl;
684 double totalvolume=0;
685 for (std::map<Material*,double>::iterator it = volume.begin();
686 it != volume.end(); ++it)
687 totalvolume += it->second;
690 if (satnum.empty()) {
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);
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;
704 if (upscaledRho > 0) {
705 upscaledRho /= totalvolume;
706 std::cout <<
"Upscaled density: " << upscaledRho << std::endl;
710IMPL_FUNC(
void, loadMaterialsFromRocklist(
const std::string& file,
711 const std::string& rocklist))
713 std::vector< std::shared_ptr<Material> > cache;
716 f.open(rocklist.c_str());
719 for (
int i=0;i<mats;++i) {
728 for (
size_t i=0;i<cache.size();++i)
729 Emin = std::min(Emin,
static_cast<Isotropic*
>(cache[i].get())->getE());
730 for (
size_t i=0;i<cache.size();++i) {
731 double lE =
static_cast<Isotropic*
>(cache[i].get())->getE();
732 static_cast<Isotropic*
>(cache[i].get())->setE(lE*Escale/Emin);
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]);
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);
750 auto deck = parser.parseFile(file);
751 std::vector<int> satnum = deck[
"SATNUM"].back().getIntData();
752 std::vector<int> cells = gv.globalCell();
753 for (
size_t i=0;i<cells.size();++i) {
755 if (satnum[k]-1 >= (
int)cache.size()) {
756 std::cerr <<
"Material " << satnum[k] <<
" referenced but not available. Check your rocklist." << std::endl;
759 materials.push_back(cache[satnum[k]-1]);
760 volume[satnum[k]-1] += gv.cellVolume(i);
763 std::cout <<
"Number of materials: " << cache.size() << std::endl;
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);
773IMPL_FUNC(
void, fixCorners(
const double* min,
const double* max))
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];
790IMPL_FUNC(
void, periodicBCs(
const double* min,
const double* max))
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;
811 periodicPlane(X,XYZ,slave[0],master[0]);
812 periodicPlane(Y,XYZ,slave[1],master[1]);
813 periodicPlane(Z,XYZ,slave[2],master[2]);
816IMPL_FUNC(
void, periodicBCsMortar(
const double* min,
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;
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));
856 std::cout <<
"\testablished YZ multiplier grid with " << n2 <<
"x1" <<
" elements" << std::endl;
860 fmin[0] = min[1]; fmin[1] = min[2];
861 fmax[0] = max[1]; fmax[1] = max[2];
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;
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);
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);
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);
891 template<
class M,
class A>
892static void applyMortarBlock(
int i,
const Matrix& B, M& T,
900 MortarBlockEvaluator<Dune::Preconditioner<Vector,Vector> > pre(upre, B);
903 for (
size_t j=0; j < B.M(); ++j)
907IMPL_FUNC(
void, setupSolvers(
const LinSolParams& params))
909 int siz = A.getOperator().N();
912 numsolvers = omp_get_max_threads();
915 if (params.type == ITERATIVE) {
916 op.reset(
new Operator(A.getOperator()));
918 upre.push_back(PC::setup(params.steps[0], params.steps[1],
919 params.coarsen_target, params.zcells,
927 if (params.mortarpre >= SCHUR) {
928 Dune::DynamicMatrix<double> T(B.M(), B.M());
929 std::cout <<
"\tBuilding preconditioner for multipliers..." << std::endl;
932 std::vector< std::shared_ptr<Dune::Preconditioner<Vector,Vector> > > pc;
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)) {
944 for (
size_t i=1;i<B.M();++i)
945 pc[i].reset(
new PCType(*upre[0]));
947 }
else if (params.mortarpre == SCHURAMG) {
948 std::shared_ptr<AMG1<SSORSmoother>::type> mpc;
951 params.coarsen_target,
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;
960 params.coarsen_target,
963 }
else if (params.mortarpre == SCHURTWOLEVEL) {
964 std::shared_ptr<AMG2Level<SSORSmoother>::type> mpc;
967 params.coarsen_target,
970 for (
size_t i=1;i<B.M();++i)
971 pc[i].reset(
new AMG2Level<SSORSmoother>::type(*mpc));
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]);
978 help.log(help.group(t).second,
979 "\t\t... still processing ... multiplier ");
982 }
else if (params.mortarpre == SIMPLE) {
987 for (Matrix::ConstRowIterator it = A.getOperator().begin();
988 it != A.getOperator().end(); ++it, ++row) {
990 for (Matrix::ConstColIterator it2 = it->begin();
991 it2 != it->end(); ++it2) {
992 if (it2.index() != row)
995 D[row][row] = 1.0/(A.getOperator()[row][row]/alpha);
1000 Dune::matMultMat(t1, D, B);
1002 Dune::transposeMatMultMat(P, B, t1);
1006 std::shared_ptr<Dune::InverseOperator<Vector,Vector> > innersolver;
1008 innersolver.reset(
new Dune::CGSolver<Vector>(*op, *upre[0], params.tol,
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,
1017 verbose?2:(params.report?1:0)));
1018 tsolver.push_back(SolverPtr(
new UzawaSolver<Vector, Vector>(innersolver, outersolver, B)));
1020 for (
int i=0;i<numsolvers;++i) {
1022 upre.push_back(PCPtr(
new PCType(*upre[0])));
1023 tmpre.push_back(MortarAmgPtr(
new MortarSchurPre<PCType>(P, B,
1025 params.symmetric)));
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],
1033 verbose?2:(params.report?1:0))));
1035 for (
int i=0;i<numsolvers;++i)
1036 tsolver.push_back(SolverPtr(
new Dune::RestartedGMResSolver<Vector>(*meval, *tmpre[i],
1040 verbose?2:(params.report?1:0))));
1044 for (
int i=0;i<numsolvers;++i) {
1046 upre.push_back(PCPtr(
new PCType(*upre[0])));
1047 tsolver.push_back(SolverPtr(
new Dune::CGSolver<Vector>(*op, *upre[copy?i:0],
1050 verbose?2:(params.report?1:0))));
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());
1063 std::cerr <<
"No direct solver available" << std::endl;
1066 siz = A.getOperator().N();
1069 for (
int i=0;i<6;++i)
1073IMPL_FUNC(
void, solve(
int loadcase))
1076 Dune::InverseOperatorResult r;
1077 u[loadcase].resize(b[loadcase].size());
1081 solver = omp_get_thread_num();
1084 tsolver[solver]->apply(u[loadcase], b[loadcase], r);
1086 std::cout <<
"\tsolution norm: " << u[loadcase].two_norm() << std::endl;
1087 }
catch (
const Dune::ISTLError& e) {
1088 std::cerr <<
"exception thrown " << e << std::endl;
BoundaryGrid::Vertex maxXminY(const std::vector< BoundaryGrid::Vertex > &in)
Find the vertex in the vector with maximum X and minimum Y.
Definition boundarygrid.cpp:361
BoundaryGrid::Vertex minXminY(const std::vector< BoundaryGrid::Vertex > &in)
Find the vertex in the vector with minimum X and minimum Y.
Definition boundarygrid.cpp:351
BoundaryGrid::Vertex maxXmaxY(const std::vector< BoundaryGrid::Vertex > &in)
Find the vertex in the vector with maximum X and maximum Y.
Definition boundarygrid.cpp:371
BoundaryGrid::Vertex minXmaxY(const 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 > ¶ms)
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:49
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 ©)
Setup preconditioner.
Definition elasticity_preconditioners.hpp:117
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 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