12 #ifndef OPM_ELASTICITY_UPSCALE_IMPL_HPP
13 #define OPM_ELASTICITY_UPSCALE_IMPL_HPP
22 namespace Elasticity {
25 #define IMPL_FUNC(A,B) template<class GridType, class PC> \
26 A ElasticityUpscale<GridType, PC>::B
31 std::vector<BoundaryGrid::Vertex> result;
32 const LeafVertexIterator
itend = gv.leafGridView().template end<dim>();
35 Dune::LeafMultipleCodimMultipleGeomTypeMapper<GridType,
36 Dune::MCMGVertexLayout> mapper(gv);
38 LeafVertexIterator
start = gv.leafGridView().template begin<dim>();
39 for (LeafVertexIterator it = start; it !=
itend; ++it) {
40 if (isOnPlane(dir,it->geometry().corner(0),coord)) {
42 v.
i = mapper.map(*it);
57 static const int V1[3][4] = {{0,2,4,6},
60 static const int V2[3][4] = {{1,3,5,7},
63 const LeafIndexSet& set = gv.leafGridView().indexSet();
70 std::map<double, std::vector<BoundaryGrid::Quad> > nodeMap;
71 for (LeafIterator cell = gv.leafGridView().template begin<0>();
72 cell != gv.leafGridView().template end<0>(); ++cell, ++c) {
73 std::vector<BoundaryGrid::Vertex> verts;
76 idx = set.subIndex(*cell,V1[i][0],dim);
77 else if (side == RIGHT)
78 idx = set.subIndex(*cell,V2[i][0],dim);
79 Dune::FieldVector<double, 3> pos = gv.vertexPosition(idx);
80 if (isOnPlane(dir,pos,coord)) {
81 for (
int j=0;
j<4;++
j) {
83 idx = set.subIndex(*cell,V1[i][
j],dim);
85 idx = set.subIndex(*cell,V2[i][j],dim);
86 pos = gv.vertexPosition(idx);
87 if (!isOnPlane(dir,pos,coord))
95 if (verts.size() == 4) {
106 std::map<double, std::vector<BoundaryGrid::Quad> >::iterator it;
107 for (it = nodeMap.begin(); it != nodeMap.end(); ++it) {
108 if (fabs(it->first-q.
v[0].
c[0]) < 1.e-7) {
109 it->second.push_back(q);
113 if (it == nodeMap.end())
114 nodeMap[q.
v[0].
c[0]].push_back(q);
121 std::map<double, std::vector<BoundaryGrid::Quad> >::const_iterator it;
122 for (it = nodeMap.begin(); it != nodeMap.end(); ++it, ++p) {
123 for (
size_t i=0;i<it->second.size();++i)
132 master.push_back(extractMasterFace(
X,min[0]));
133 master.push_back(extractMasterFace(
Y,min[1]));
134 master.push_back(extractMasterFace(
Z,min[2]));
136 slave.push_back(extractFace(
X,max[0]));
137 slave.push_back(extractFace(
Y,max[1]));
138 slave.push_back(extractFace(
Z,max[2]));
141 IMPL_FUNC(
void, findBoundaries(
double* min,
double* max))
143 max[0] = max[1] = max[2] = -1e5;
144 min[0] = min[1] = min[2] = 1e5;
145 const LeafVertexIterator
itend = gv.leafGridView().template end<dim>();
148 LeafVertexIterator
start = gv.leafGridView().template begin<dim>();
149 for (LeafVertexIterator it = start; it !=
itend; ++it) {
150 for (
int i=0;i<3;++i) {
151 min[i] =
std::min(min[i],it->geometry().corner(0)[i]);
152 max[i] = std::max(max[i],it->geometry().corner(0)[i]);
160 MPC* mpc =
new MPC(slave,log2(dir)+1);
165 for (
int i=0;i<4;++i)
173 const std::vector<BoundaryGrid::Vertex>& slave,
176 for (
size_t i=0;i<slave.size();++i) {
178 if (master.
find(coord,slave[i])) {
179 addMPC(
X,slave[i].i,coord);
180 addMPC(
Y,slave[i].i,coord);
181 addMPC(
Z,slave[i].i,coord);
195 std::vector<std::vector<int> > nodes;
196 nodes.resize(b.
size());
202 nodes[0].push_back(-1);
203 for (
size_t e=0; e < b.
size(); ++e) {
205 for (
int i2=0; i2 < n2; ++i2) {
207 nodes[e].push_back(nodes[e-1][i2*n1+n1-1]);
209 int start = (e==0 && i2 != 0)?0:1;
212 if (i2 == n2-1 && n2 > 2) {
213 for (
int i1=(e==0?0:1); i1 < n1; ++i1) {
214 nodes[e].push_back(nodes[e][i1]);
217 for (
int i1=start; i1 < n1; ++i1) {
221 nodes[e].push_back(totalDOFs++);
232 int dir,
int n1,
int n2,
237 std::vector<std::vector<int> > lnodes =
renumber(interface, n1+1,
240 Bpatt.resize(A.getEqns());
243 for (
size_t p=0;p<interface.
size();++p) {
244 for (
size_t q=0;q<b1.
colSize(p);++q) {
245 for (
size_t i=0;i<4;++i) {
246 for (
size_t d=0;d<3;++d) {
252 for (
size_t j=0;
j<lnodes[p].size();++
j) {
253 int indexj = 3*lnodes[p][
j]+d;
255 Bpatt[dof].insert(indexj+colofs);
260 int dof = A.getEquationForDof(b1.
getQuad(p,q).
v[i].
i,d);
262 for (
size_t j=0;
j<lnodes[p].size();++
j) {
263 int indexj = 3*lnodes[p][
j]+d;
265 Bpatt[dof].insert(indexj+colofs);
292 std::vector<std::vector<int> > lnodes =
renumber(interface, n1+1,
295 Dune::GeometryType gt;
298 int quadorder = std::max((1.0+n1+0.5)/2.0,(1.0+n2+0.5)/2.0);
299 quadorder = std::max(quadorder, 2);
300 const Dune::QuadratureRule<ctype,2>& rule =
301 Dune::QuadratureRules<ctype,2>::rule(gt,quadorder);
304 typename Dune::QuadratureRule<ctype,2>::const_iterator r;
305 Dune::DynamicMatrix<ctype> E(ubasis.
n,(n1+1)*(n2+1),0.0);
307 for (
int g=0;g<5;++g) {
308 for (
int p=help.group(g).first;p<help.group(g).second;++p) {
311 for (
size_t q=0;q<b1.
colSize(p);++q) {
315 for (r = rule.begin(); r != rule.end();++r) {
316 ctype detJ = hex.integrationElement(r->position());
321 lg.local(hex.global(r->position()));
322 assert(loc[0] <= 1.0+1.e-4 && loc[0] >= 0.0 && loc[1] <= 1.0+1.e-4 && loc[1] >= 0.0);
323 for (
int i=0;i<ubasis.
n;++i) {
324 for (
int j=0;
j<lbasis.
size();++
j) {
325 E[i][
j] += ubasis[i].evaluateFunction(r->position())*
326 lbasis[
j].evaluateFunction(loc)*detJ*r->weight();
332 for (
int d=0;d<3;++d) {
333 for (
int i=0;i<4;++i) {
334 const MPC* mpc = A.getMPC(qu.
v[i].
i,d);
339 for (
size_t j=0;
j<lnodes[p].size();++
j) {
340 int indexj = lnodes[p][
j]*3+d;
342 B[indexi][indexj+colofs] += alpha*E[i][
j];
347 int indexi = A.getEquationForDof(qu.
v[i].
i,d);
349 for (
size_t j=0;
j<lnodes[p].size();++
j) {
350 int indexj = lnodes[p][
j]*3+d;
352 B[indexi][indexj+colofs] += alpha*E[i][
j];
360 help.log(g,
"\t\t\t... still processing ... pillar ");
365 GlobalCoordinate coord,
366 const NodeValue& value))
368 typedef typename GridType::LeafGridView::template Codim<dim>::Iterator VertexLeafIterator;
369 const VertexLeafIterator itend = gv.leafGridView().template end<dim>();
372 Dune::LeafMultipleCodimMultipleGeomTypeMapper<GridType,
373 Dune::MCMGVertexLayout> mapper(gv);
376 for (VertexLeafIterator it = gv.leafGridView().template begin<dim>(); it !=
itend; ++it) {
377 if (isOnPoint(it->geometry().corner(0),coord)) {
378 int indexi = mapper.map(*it);
379 A.updateFixedNode(indexi,std::make_pair(dir,value));
385 GlobalCoordinate coord,
388 if (plane < X || plane >
Z)
391 ctype delta = fabs(value-coord[p]);
397 const NodeValue& value))
399 typedef typename GridType::LeafGridView::template Codim<dim>::Iterator VertexLeafIterator;
400 const VertexLeafIterator itend = gv.leafGridView().template end<dim>();
403 Dune::LeafMultipleCodimMultipleGeomTypeMapper<GridType,
404 Dune::MCMGVertexLayout> mapper(gv);
407 for (VertexLeafIterator it = gv.leafGridView().template begin<dim>(); it !=
itend; ++it) {
408 if (isOnLine(dir,it->geometry().corner(0),x,y)) {
409 int indexi = mapper.map(*it);
410 A.updateFixedNode(indexi,std::make_pair(
XYZ,value));
416 GlobalCoordinate coord,
419 if (dir < X || dir >
Z)
421 int ix = int(log2(dir)+1) % 3;
422 int iy = int(log2(dir)+2) % 3;
423 ctype delta = x-coord[ix];
424 if (delta > tol || delta < -tol)
427 if (delta > tol || delta < -tol)
434 GlobalCoordinate point))
436 GlobalCoordinate delta = point-coord;
437 return delta.one_norm() < tol;
440 IMPL_FUNC(
void, assemble(
int loadcase,
bool matrix))
442 const int comp = 3+(dim-2)*3;
443 static const int bfunc = 4+(dim-2)*4;
445 Dune::FieldVector<ctype,comp>
eps0;
454 for (
int i=0;i<2;++i) {
455 if (color[1].size() && matrix)
456 std::cout <<
"\tprocessing " << (i==0?
"red ":
"black ") <<
"elements" << std::endl;
457 #pragma omp parallel for schedule(static)
458 for (
size_t j=0;
j<color[i].size();++
j) {
459 Dune::FieldMatrix<ctype,comp,comp> C;
460 Dune::FieldMatrix<ctype,dim*bfunc,dim*bfunc> K;
461 Dune::FieldMatrix<ctype,dim*bfunc,dim*bfunc>* KP=0;
462 Dune::FieldVector<ctype,dim*bfunc> ES;
463 Dune::FieldVector<ctype,dim*bfunc>* EP=0;
469 for (
size_t k=0;k<color[i][
j].size();++k) {
470 LeafIterator it = gv.leafGridView().template begin<0>();
471 for (
int l=0;l<color[i][
j][k];++l)
473 materials[color[i][
j][k]]->getConstitutiveMatrix(C);
475 Dune::GeometryType gt = it->type();
477 Dune::FieldMatrix<ctype,dim*bfunc,dim*bfunc> Aq;
482 const Dune::QuadratureRule<ctype,dim>& rule = Dune::QuadratureRules<ctype,dim>::rule(gt,2);
483 for (
typename Dune::QuadratureRule<ctype,dim>::const_iterator r = rule.begin();
484 r != rule.end() ; ++r) {
486 Dune::FieldMatrix<ctype,dim,dim> jacInvTra =
487 it->geometry().jacobianInverseTransposed(r->position());
489 ctype detJ = it->geometry().integrationElement(r->position());
490 if (detJ <= 1.e-5 && verbose) {
491 std::cout <<
"cell " << color[i][
j][k] <<
" is (close to) degenerated, detJ " << detJ << std::endl;
493 for (
int i=0;i<4;++i)
494 zdiff = std::max(zdiff, it->geometry().corner(i+4)[2]-it->geometry().corner(i)[2]);
495 std::cout <<
" - Consider setting ctol larger than " << zdiff << std::endl;
498 Dune::FieldMatrix<ctype,comp,dim*bfunc> B;
499 E.getBmatrix(B,r->position(),jacInvTra);
502 E.getStiffnessMatrix(Aq,B,C,detJ*r->weight());
508 Dune::FieldVector<ctype,dim*bfunc> temp;
509 temp = Dune::FMatrixHelp::multTransposed(B,Dune::FMatrixHelp::mult(C,eps0));
510 temp *= -detJ*r->weight();
514 A.addElement(KP,EP,it,(loadcase > -1)?&b[loadcase]:NULL);
521 averageStress(Dune::FieldVector<ctype,comp>& sigma,
522 const Vector& u,
int loadcase))
524 if (loadcase < 0 || loadcase > 5)
527 static const int bfunc = 4+(dim-2)*4;
529 const LeafIterator itend = gv.leafGridView().template end<0>();
531 Dune::FieldMatrix<ctype,comp,comp> C;
532 Dune::FieldVector<ctype,comp>
eps0;
538 for (LeafIterator it = gv.leafGridView().template begin<0>(); it !=
itend; ++it) {
539 materials[m++]->getConstitutiveMatrix(C);
541 Dune::GeometryType gt = it->type();
543 Dune::FieldVector<ctype,bfunc*dim> v;
544 A.extractValues(v,u,it);
547 const Dune::QuadratureRule<ctype,dim>& rule = Dune::QuadratureRules<ctype,dim>::rule(gt,2);
548 for (
typename Dune::QuadratureRule<ctype,dim>::const_iterator r = rule.begin();
549 r != rule.end() ; ++r) {
551 Dune::FieldMatrix<ctype,dim,dim> jacInvTra =
552 it->geometry().jacobianInverseTransposed(r->position());
554 ctype detJ = it->geometry().integrationElement(r->position());
556 volume += detJ*r->weight();
558 Dune::FieldMatrix<ctype,comp,dim*bfunc> B;
559 E.getBmatrix(B,r->position(),jacInvTra);
561 Dune::FieldVector<ctype,comp> s;
562 E.getStressVector(s,v,eps0,B,C);
563 s *= detJ*r->weight();
569 sigma /= Escale/Emin;
572 IMPL_FUNC(
void, loadMaterialsFromGrid(
const std::string& file))
574 typedef std::map<std::pair<double,double>,
575 std::shared_ptr<Material> > MaterialMap;
583 if (file ==
"uniform") {
584 int cells = gv.size(0);
585 Emod.insert(Emod.begin(),
cells,100.f);
586 Poiss.insert(Poiss.begin(),
cells,0.38f);
588 Opm::ParserPtr parser(
new Opm::Parser());
590 Opm::DeckConstPtr
deck(parser->parseFile(file , parseMode));
591 if (
deck->hasKeyword(
"YOUNGMOD") &&
deck->hasKeyword(
"POISSONMOD")) {
592 Emod =
deck->getKeyword(
"YOUNGMOD")->getRawDoubleData();
593 Poiss =
deck->getKeyword(
"POISSONMOD")->getRawDoubleData();
594 std::vector<double>::const_iterator it = std::min_element(Poiss.begin(), Poiss.end());
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->getKeyword(
"LAMEMOD")->getRawDoubleData();
602 std::vector<double> shear =
deck->getKeyword(
"SHEARMOD")->getRawDoubleData();
603 Emod.resize(lame.size());
604 Poiss.resize(lame.size());
605 for (
size_t i=0;i<lame.size();++i) {
606 Emod[i] = shear[i]*(3*lame[i]+2*shear[i])/(lame[i]+shear[i]);
607 Poiss[i] = 0.5*lame[i]/(lame[i]+shear[i]);
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->getKeyword(
"BULKMOD")->getRawDoubleData();
618 std::vector<double> shear =
deck->getKeyword(
"SHEARMOD")->getRawDoubleData();
619 Emod.resize(bulk.size());
620 Poiss.resize(bulk.size());
621 for (
size_t i=0;i<bulk.size();++i) {
622 Emod[i] = 9*bulk[i]*shear[i]/(3*bulk[i]+shear[i]);
623 Poiss[i] = 0.5*(3*bulk[i]-2*shear[i])/(3*bulk[i]+shear[i]);
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->getKeyword(
"PERMX")->getRawDoubleData();
635 Poiss =
deck->getKeyword(
"PORO")->getRawDoubleData();
637 std::cerr <<
"No material data found in eclipse file, aborting" << std::endl;
640 if (
deck->hasKeyword(
"SATNUM"))
641 satnum =
deck->getKeyword(
"SATNUM")->getIntData();
642 if (
deck->hasKeyword(
"RHO"))
643 rho =
deck->getKeyword(
"RHO")->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();
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);
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 j=0; j < volumeFractions.size(); ++
j) {
699 volumeFractions[
j] /= totalvolume;
700 std::cout <<
"SATNUM " << j+1 <<
": " << volumeFractions[
j]*100 <<
'%' << std::endl;
706 std::cout <<
"Upscaled density: " <<
upscaledRho << std::endl;
710 IMPL_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)
730 for (
size_t i=0;i<cache.size();++i) {
731 double E = ((
Isotropic*)cache[i].
get())->getE();
732 ((
Isotropic*)cache[i].
get())->setE(E*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);
749 Opm::ParserPtr parser(
new Opm::Parser());
751 Opm::DeckConstPtr
deck(parser->parseFile(file , parseMode));
752 std::vector<int> satnum =
deck->getKeyword(
"SATNUM")->getIntData();
753 std::vector<int> cells = gv.globalCell();
754 for (
size_t i=0;i<cells.size();++i) {
756 if (satnum[k]-1 >= (
int)cache.size()) {
757 std::cerr <<
"Material " << satnum[k] <<
" referenced but not available. Check your rocklist." << std::endl;
760 materials.push_back(cache[satnum[k]-1]);
761 volume[satnum[k]-1] += gv.cellVolume(i);
764 std::cout <<
"Number of materials: " << cache.size() << std::endl;
766 double totalvolume = std::accumulate(volume.begin(),volume.end(),0.f);
767 for (
size_t i=0;i<cache.size();++i) {
768 std::cout <<
" SATNUM " << i+1 <<
": " << 100.f*volume[i]/totalvolume <<
'%' << std::endl;
769 volumeFractions.push_back(volume[i]/totalvolume);
774 IMPL_FUNC(
void, fixCorners(
const double* min,
const double* max))
776 ctype c[8][3] = {{min[0],min[1],min[2]},
777 {max[0],min[1],min[2]},
778 {min[0],max[1],min[2]},
779 {max[0],max[1],min[2]},
780 {min[0],min[1],max[2]},
781 {max[0],min[1],max[2]},
782 {min[0],max[1],max[2]},
783 {max[0],max[1],max[2]}};
784 for (
int i=0;i<8;++i) {
785 GlobalCoordinate coord;
786 coord[0] = c[i][0]; coord[1] = c[i][1]; coord[2] = c[i][2];
791 IMPL_FUNC(
void, periodicBCs(
const double* min,
const double* max))
804 std::cout <<
"Xslave " << slave[0].size() <<
" "
805 <<
"Yslave " << slave[1].size() <<
" "
806 <<
"Zslave " << slave[2].size() << std::endl;
807 std::cout <<
"Xmaster " << master[0].size() <<
" "
808 <<
"Ymaster " << master[1].size() <<
" "
809 <<
"Zmaster " << master[2].size() << std::endl;
812 periodicPlane(
X,
XYZ,slave[0],master[0]);
813 periodicPlane(
Y,
XYZ,slave[1],master[1]);
814 periodicPlane(
Z,
XYZ,slave[2],master[2]);
817 IMPL_FUNC(
void, periodicBCsMortar(
const double* min,
838 std::cout <<
"\textracting nodes on top face..." << std::endl;
839 slave.push_back(extractFace(
Z,max[2]));
840 std::cout <<
"\treconstructing bottom face..." << std::endl;
841 BoundaryGrid bottom = extractMasterFace(
Z,min[2]);
842 std::cout <<
"\testablishing couplings on top/bottom..." << std::endl;
843 periodicPlane(
Z,
XYZ,slave[0],bottom);
844 std::cout <<
"\tinitializing matrix..." << std::endl;
848 std::cout <<
"\treconstructing left face..." << std::endl;
849 master.push_back(extractMasterFace(
X, min[0], LEFT,
true));
850 std::cout <<
"\treconstructing right face..." << std::endl;
851 master.push_back(extractMasterFace(
X, max[0], RIGHT,
true));
852 std::cout <<
"\treconstructing front face..." << std::endl;
853 master.push_back(extractMasterFace(
Y, min[1], LEFT,
true));
854 std::cout <<
"\treconstructing back face..." << std::endl;
855 master.push_back(extractMasterFace(
Y, max[1], RIGHT,
true));
857 std::cout <<
"\testablished YZ multiplier grid with " << n2 <<
"x1" <<
" elements" << std::endl;
861 fmin[0] = min[1]; fmin[1] = min[2];
862 fmax[0] = max[1]; fmax[1] = max[2];
865 fmin[0] = min[0]; fmin[1] = min[2];
866 fmax[0] = max[0]; fmax[1] = max[2];
867 std::cout <<
"\testablished XZ multiplier grid with " << n1 <<
"x1" <<
" elements" << std::endl;
870 addBBlockMortar(master[0], lambdax, 0, 1, p2, 0);
871 int eqns = addBBlockMortar(master[1], lambdax, 0, 1, p2, 0);
872 addBBlockMortar(master[2], lambday, 1, 1, p1, eqns);
873 int eqns2 = addBBlockMortar(master[3], lambday, 1, 1, p1, eqns);
879 std::cout <<
"\tassembling YZ mortar matrix..." << std::endl;
880 assembleBBlockMortar(master[0], lambdax, 0, 1, p2, 0);
881 assembleBBlockMortar(master[1], lambdax, 0, 1, p2, 0, -1.0);
884 std::cout <<
"\tassembling XZ mortar matrix..." << std::endl;
885 assembleBBlockMortar(master[2], lambday, 1, 1, p1, eqns);
886 assembleBBlockMortar(master[3], lambday, 1, 1, p1, eqns, -1.0);
892 template<
class M,
class A>
904 for (
size_t j=0; j < B.M(); ++
j)
908 IMPL_FUNC(
void, setupSolvers(
const LinSolParams& params))
910 int siz = A.getOperator().N();
913 numsolvers = omp_get_max_threads();
917 op.reset(
new Operator(A.getOperator()));
919 upre.push_back(PC::setup(params.steps[0], params.steps[1],
920 params.coarsen_target, params.zcells,
928 if (params.mortarpre >=
SCHUR) {
929 Dune::DynamicMatrix<double> T(B.M(), B.M());
930 std::cout <<
"\tBuilding preconditioner for multipliers..." << std::endl;
933 std::vector< std::shared_ptr<Dune::Preconditioner<Vector,Vector> > > pc;
936 if (params.mortarpre ==
SCHUR ||
938 params.pre ==
AMG) ||
945 for (
size_t i=1;i<B.M();++i)
946 pc[i].reset(
new PCType(*upre[0]));
948 }
else if (params.mortarpre ==
SCHURAMG) {
949 std::shared_ptr<AMG1<SSORSmoother>::type> mpc;
950 pc[0] = mpc = AMG1<SSORSmoother>::setup(params.steps[0],
952 params.coarsen_target,
955 for (
size_t i=1;i<B.M();++i)
956 pc[i].reset(
new AMG1<SSORSmoother>::type(*mpc));
958 std::shared_ptr<Schwarz::type> mpc;
959 pc[0] = mpc = Schwarz::setup(params.steps[0],
961 params.coarsen_target,
965 std::shared_ptr<AMG2Level<SSORSmoother>::type> mpc;
966 pc[0] = mpc = AMG2Level<SSORSmoother>::setup(params.steps[0],
968 params.coarsen_target,
971 for (
size_t i=1;i<B.M();++i)
972 pc[i].reset(
new AMG2Level<SSORSmoother>::type(*mpc));
974 for (
int t=0;t<5;++t) {
975 #pragma omp parallel for schedule(static)
976 for (
int i=help.group(t).first; i < help.group(t).second; ++i)
979 help.log(help.group(t).second,
980 "\t\t... still processing ... multiplier ");
983 }
else if (params.mortarpre ==
SIMPLE) {
988 for (Matrix::ConstRowIterator it = A.getOperator().begin();
989 it != A.getOperator().end(); ++it, ++row) {
991 for (Matrix::ConstColIterator it2 = it->begin();
992 it2 != it->end(); ++it2) {
993 if (it2.index() != row)
996 D[row][row] = 1.0/(A.getOperator()[row][row]/alpha);
1001 Dune::matMultMat(t1, D, B);
1003 Dune::transposeMatMultMat(P, B, t1);
1007 std::shared_ptr<Dune::InverseOperator<Vector,Vector> > innersolver;
1009 innersolver.reset(
new Dune::CGSolver<Vector>(*op, *upre[0], params.tol,
1011 verbose?2:(params.report?1:0)));
1012 op2.reset(
new SchurEvaluator(*innersolver, B));
1013 lprep.reset(
new LUSolver(P));
1014 lpre.reset(
new SeqLU(*lprep));
1015 std::shared_ptr<Dune::InverseOperator<Vector,Vector> > outersolver;
1016 outersolver.reset(
new Dune::CGSolver<Vector>(*op2, *lpre, params.tol*10,
1018 verbose?2:(params.report?1:0)));
1023 upre.push_back(PCPtr(
new PCType(*upre[0])));
1026 params.symmetric)));
1029 if (params.symmetric) {
1031 tsolver.push_back(SolverPtr(
new Dune::MINRESSolver<Vector>(*meval, *tmpre[i],
1034 verbose?2:(params.report?1:0))));
1037 tsolver.push_back(SolverPtr(
new Dune::RestartedGMResSolver<Vector>(*meval, *tmpre[i],
1041 verbose?2:(params.report?1:0),
true)));
1047 upre.push_back(PCPtr(
new PCType(*upre[0])));
1048 tsolver.push_back(SolverPtr(
new Dune::CGSolver<Vector>(*op, *upre[copy?i:0],
1051 verbose?2:(params.report?1:0))));
1057 0, A.getOperator().M(),
true);
1058 #if HAVE_UMFPACK || HAVE_SUPERLU
1059 tsolver.push_back(SolverPtr(
new LUSolver(A.getOperator(),
1060 verbose?2:(params.report?1:0))));
1062 tsolver.push_back(tsolver.front());
1064 std::cerr <<
"No direct solver available" << std::endl;
1067 siz = A.getOperator().N();
1070 for (
int i=0;i<6;++i)
1077 Dune::InverseOperatorResult r;
1078 u[loadcase].resize(b[loadcase].size(),
false);
1082 solver = omp_get_thread_num();
1085 tsolver[solver]->apply(u[loadcase], b[loadcase], r);
1087 std::cout <<
"\tsolution norm: " << u[loadcase].two_norm() << std::endl;
1088 }
catch (Dune::ISTLError& e) {
1089 std::cerr <<
"exception thrown " << e << std::endl;
LeafVertexIterator start
Definition: elasticity_upscale_impl.hpp:148
diagonal approximation of A
Definition: elasticity_upscale.hpp:62
A class describing a linear, quadrilateral element.
Definition: boundarygrid.hh:81
static void extract(FaceCoord &res, const GlobalCoordinate &coord, int dir)
Helper function for extracting given 2D coordinates from a 3D coordinate.
std::vector< int > satnum
Definition: elasticity_upscale_impl.hpp:579
upscaledRho
Definition: elasticity_upscale_impl.hpp:581
std::vector< double > Emod
Definition: elasticity_upscale_impl.hpp:577
std::map< Material *, double > volume
Definition: elasticity_upscale_impl.hpp:659
schur + primary preconditioner
Definition: elasticity_upscale.hpp:63
static Matrix diagonal(size_t N)
Returns a diagonal matrix.
FaceCoord c
Coordinates of the vertex (2D)
Definition: boundarygrid.hh:66
Definition: elasticity_upscale.hpp:56
std::vector< int > cells
Definition: elasticity_upscale_impl.hpp:657
Isotropic linear elastic material.
Definition: materials.hh:25
Definition: applier.hpp:18
static Matrix fromDense(const Dune::DynamicMatrix< double > &T)
Create a sparse matrix from a dense matrix.
std::vector< double > evalBasis(double xi, double eta) const
Evaluate the basis functions.
BoundaryGrid::Vertex minXmaxY(std::vector< BoundaryGrid::Vertex > &in)
Find the vertex in the vector with minimum X and maximum Y.
int node
Node number identifying this DOF.
Definition: mpc.hh:83
MaterialMap cache
Definition: elasticity_upscale_impl.hpp:573
void addToColumn(size_t col, const Quad &quad)
Definition: boundarygrid.hh:118
int size()
Definition: shapefunctions.hpp:322
const Quad & getQuad(int col, int index) const
Definition: boundarygrid.hh:141
if(loadcase >-1)
Definition: elasticity_upscale_impl.hpp:447
Vertex v[4]
Vertices.
Definition: boundarygrid.hh:101
static Material * create(int ID, const Dune::DynamicVector< double > ¶ms)
Creates a material object of a given type.
Definition: mortar_schur_precond.hpp:37
schur + twolevel
Definition: elasticity_upscale.hpp:66
void add(const Quad &quad)
Add a quad to the grid.
std::vector< double > rho
Definition: elasticity_upscale_impl.hpp:580
std::vector< double > Poiss
Definition: elasticity_upscale_impl.hpp:578
schur + amg
Definition: elasticity_upscale.hpp:64
Definition: mortar_evaluator.hpp:27
A class for representing a general multi-point constraint equation.
Definition: mpc.hh:58
bool find(Vertex &res, const Vertex &coord) const
Locate the position of a vertex on the grid.
void apply(const Vector &x, Vector &y) const
Apply the multiplier block.
Definition: mortar_schur.hpp:47
static std::vector< std::vector< int > > renumber(const BoundaryGrid &b, int n1, int n2, int &totalDOFs)
Static helper to renumber a Q4 mesh to a P_1/P_N lag mesh.
Definition: elasticity_upscale_impl.hpp:191
static const int bfunc
Definition: elasticity_upscale_impl.hpp:443
Definition: mortar_schur.hpp:27
const DOF & getMaster(size_t i) const
Returns a reference to the i'th master DOF.
Definition: mpc.hh:141
int i
Index of the vertex.
Definition: boundarygrid.hh:63
A class describing a 2D vertex.
Definition: boundarygrid.hh:57
Dune::BCRSMatrix< Dune::FieldMatrix< double, 1, 1 > > Matrix
A sparse matrix holding our operator.
Definition: matrixops.hpp:23
size_t getNoMaster() const
Returns the number of master DOFs.
Definition: mpc.hh:144
min[0]
Definition: elasticity_upscale_impl.hpp:144
static void fromAdjacency(Matrix &A, const AdjacencyPattern &adj, int rows, int cols)
Create a sparse matrix from a given adjacency pattern.
Geometry class for general hexagons.
Definition: boundarygrid.hh:283
static const P1ShapeFunctionSet & instance()
Get the only instance of this class.
Definition: shapefunctions.hpp:193
size_t size() const
Obtain the number of quads in the grid.
Definition: boundarygrid.hh:147
A class describing a quad grid.
Definition: boundarygrid.hh:28
BoundaryGrid::Vertex maxXminY(std::vector< BoundaryGrid::Vertex > &in)
Find the vertex in the vector with maximum X and minimum Y.
Definition: elasticity_upscale.hpp:49
static BoundaryGrid uniform(const FaceCoord &min, const FaceCoord &max, int k1, int k2, bool dc=false)
Establish an uniform quad grid.
Direction
An enum for specification of global coordinate directions.
Definition: mpc.hh:24
Dune::FieldVector< ctype, comp > eps0
Definition: elasticity_upscale_impl.hpp:445
static std::cout<< "Xslave "<< slave[0].size()<< " "<< "Yslave "<< slave[1].size()<< " "<< "Zslave "<< slave[2].size()<< std::endl;std::cout<< "Xmaster "<< master[0].size()<< " "<< "Ymaster "<< master[1].size()<< " "<< "Zmaster "<< master[2].size()<< std::endl;periodicPlane(X, XYZ, slave[0], master[0]);periodicPlane(Y, XYZ, slave[1], master[1]);periodicPlane(Z, XYZ, slave[2], master[2]);}IMPL_FUNC(void, periodicBCsMortar(const double *min, const double *max, int n1, int n2, int p1, int p2)){fixCorners(min, max);std::cout<< "\textracting nodes on top face..."<< std::endl;slave.push_back(extractFace(Z, max[2]));std::cout<< "\treconstructing bottom face..."<< std::endl;BoundaryGrid bottom=extractMasterFace(Z, min[2]);std::cout<< "\testablishing couplings on top/bottom..."<< std::endl;periodicPlane(Z, XYZ, slave[0], bottom);std::cout<< "\tinitializing matrix..."<< std::endl;A.initForAssembly();std::cout<< "\treconstructing left face..."<< std::endl;master.push_back(extractMasterFace(X, min[0], LEFT, true));std::cout<< "\treconstructing right face..."<< std::endl;master.push_back(extractMasterFace(X, max[0], RIGHT, true));std::cout<< "\treconstructing front face..."<< std::endl;master.push_back(extractMasterFace(Y, min[1], LEFT, true));std::cout<< "\treconstructing back face..."<< std::endl;master.push_back(extractMasterFace(Y, max[1], RIGHT, true));std::cout<< "\testablished YZ multiplier grid with "<< n2<< "x1"<< " elements"<< std::endl;BoundaryGrid::FaceCoord fmin, fmax;fmin[0]=min[1];fmin[1]=min[2];fmax[0]=max[1];fmax[1]=max[2];BoundaryGrid lambdax=BoundaryGrid::uniform(fmin, fmax, n2, 1, true);fmin[0]=min[0];fmin[1]=min[2];fmax[0]=max[0];fmax[1]=max[2];std::cout<< "\testablished XZ multiplier grid with "<< n1<< "x1"<< " elements"<< std::endl;BoundaryGrid lambday=BoundaryGrid::uniform(fmin, fmax, n1, 1, true);addBBlockMortar(master[0], lambdax, 0, 1, p2, 0);int eqns=addBBlockMortar(master[1], lambdax, 0, 1, p2, 0);addBBlockMortar(master[2], lambday, 1, 1, p1, eqns);int eqns2=addBBlockMortar(master[3], lambday, 1, 1, p1, eqns);MatrixOps::fromAdjacency(B, Bpatt, A.getEqns(), eqns+eqns2);Bpatt.clear();std::cout<< "\tassembling YZ mortar matrix..."<< std::endl;assembleBBlockMortar(master[0], lambdax, 0, 1, p2, 0);assembleBBlockMortar(master[1], lambdax, 0, 1, p2, 0,-1.0);std::cout<< "\tassembling XZ mortar matrix..."<< std::endl;assembleBBlockMortar(master[2], lambday, 1, 1, p1, eqns);assembleBBlockMortar(master[3], lambday, 1, 1, p1, eqns,-1.0);master.clear();slave.clear();}template< class M, class A > void applyMortarBlock(int i, const Matrix &B, M &T, A &upre)
Definition: elasticity_upscale_impl.hpp:893
Definition: shapefunctions.hpp:272
Opm::DeckConstPtr deck(parser->parseFile(file, parseMode))
Definition: uzawa_solver.hpp:24
Opm::ParseMode parseMode
Definition: elasticity_upscale_impl.hpp:589
int numsolvers
Definition: elasticity_upscale_impl.hpp:911
size_t colSize(int i) const
Definition: boundarygrid.hh:152
static Matrix augment(const Matrix &A, const Matrix &B, size_t r0, size_t c0, bool symmetric)
Augment a matrix with another.
determineSideFaces(min, max)
IMPL_FUNC(std::vector< BoundaryGrid::Vertex >, extractFace(Direction dir, ctype coord))
Definition: elasticity_upscale_impl.hpp:28
master push_back(extractMasterFace(Y, min[1]))
BoundaryGrid::Vertex maxXmaxY(std::vector< BoundaryGrid::Vertex > &in)
Find the vertex in the vector with maximum X and maximum Y.
Dune::BlockVector< Dune::FieldVector< double, 1 > > Vector
A vector holding our RHS.
Definition: matrixops.hpp:29
Singleton handler for the set of LinearShapeFunction.
Definition: shapefunctions.hpp:180
const LeafVertexIterator itend
Definition: elasticity_upscale_impl.hpp:145
Quad * q
The quad containing the vertex (if search has been done)
Definition: boundarygrid.hh:69
Helper class for progress logging during time consuming processes.
Definition: logutils.hpp:21
Definition: elasticity_upscale.hpp:53
Definition: shapefunctions.hpp:184
int j
Definition: elasticity_upscale_impl.hpp:658
Definition: elasticity_upscale.hpp:55
void addMaster(int n, int d, double c=double(1), double tol=double(1.0e-8))
Adds a master DOF to the constraint equation.
Definition: mpc.hh:113
schur + schwarz+lu
Definition: elasticity_upscale.hpp:65
BoundaryGrid::Vertex minXminY(std::vector< BoundaryGrid::Vertex > &in)
Find the vertex in the vector with minimum X and minimum Y.
Dune::FieldVector< double, 2 > FaceCoord
A coordinate on the quad grid.
Definition: boundarygrid.hh:33