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
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)) {
44 v.
i = mapper.index(*it);
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))
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)
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;
147 const LeafVertexIterator
itend = gv.leafGridView().template end<dim>();
150 LeafVertexIterator
start = gv.leafGridView().template begin<dim>();
152 for (
int i=0;i<3;++i) {
154 max[i] = std::max(max[i],it->geometry().corner(0)[i]);
162 MPC* mpc =
new MPC(slavenode,log2(
float(dir))+1);
167 for (
int i=0;i<4;++i)
175 const std::vector<BoundaryGrid::Vertex>& slavepointgrid,
178 for (
size_t i=0;i<slavepointgrid.size();++i) {
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);
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) {
223 nodes[e].push_back(totalDOFs++);
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) {
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);
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) {
312 for (
size_t q=0;q<b1.
colSize(p);++q) {
316 for (r = rule.begin(); r != rule.end();++r) {
317 ctype detJ = hex.integrationElement(r->position());
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);
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 ");
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));
386 GlobalCoordinate coord,
389 if (plane < X || plane >
Z)
391 int p = log2(
float(plane));
392 ctype delta = fabs(value-coord[p]);
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));
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)
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);
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;
584 if (file ==
"uniform") {
585 int cells = gv.size(0);
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;
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();
660 for (
size_t i=0;i<
cells.size();++i) {
662 MaterialMap::iterator it;
664 assert(gv.cellVolume(i) > 0);
665 volume[it->second.get()] += gv.cellVolume(i);
666 materials.push_back(it->second);
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);
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();
687 totalvolume += it->second;
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;
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)
730 for (
size_t i=0;i<
cache.size();++i) {
735 std::vector<double>
volume;
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) {
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) {
756 std::cerr <<
"Material " <<
satnum[k] <<
" referenced but not available. Check your rocklist." << std::endl;
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))
778 {max[0],max[1],
min[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))
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]);
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>
903 for (
size_t j=0;
j < B.M(); ++
j)
907IMPL_FUNC(
void, setupSolvers(
const LinSolParams& params))
909 int siz = A.getOperator().N();
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 ||
937 params.pre ==
AMG) ||
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)
957 std::shared_ptr<Schwarz::type> mpc;
960 params.coarsen_target,
964 std::shared_ptr<AMG2Level<SSORSmoother>::type> mpc;
967 params.coarsen_target,
970 for (
size_t i=1;i<B.M();++i)
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)
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)));
1022 upre.push_back(PCPtr(
new PCType(*upre[0])));
1025 params.symmetric)));
1028 if (params.symmetric) {
1030 tsolver.push_back(SolverPtr(
new Dune::MINRESSolver<Vector>(*meval, *tmpre[i],
1033 verbose?2:(params.report?1:0))));
1036 tsolver.push_back(SolverPtr(
new Dune::RestartedGMResSolver<Vector>(*meval, *tmpre[i],
1040 verbose?2:(params.report?1:0))));
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))));
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)
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 (Dune::ISTLError& e) {
1088 std::cerr <<
"exception thrown " << e << std::endl;
Helper class for progress logging during time consuming processes.
Definition: logutils.hpp:21
void log(int it, const std::string &prefix)
Log to the terminal.
Definition: logutils.hpp:54
std::pair< int, int > group(int group)
Returns the start and end offsets of a chunk group.
Definition: logutils.hpp:46
A class describing a linear, quadrilateral element.
Definition: boundarygrid.hh:86
Vertex v[4]
Vertices.
Definition: boundarygrid.hh:106
std::vector< double > evalBasis(double xi, double eta) const
Evaluate the basis functions.
A class describing a 2D vertex.
Definition: boundarygrid.hh:62
Quad * q
The quad containing the vertex (if search has been done)
Definition: boundarygrid.hh:74
FaceCoord c
Coordinates of the vertex (2D)
Definition: boundarygrid.hh:71
int i
Index of the vertex.
Definition: boundarygrid.hh:68
A class describing a quad grid.
Definition: boundarygrid.hh:33
void addToColumn(size_t col, const Quad &quad)
Definition: boundarygrid.hh:123
static void extract(FaceCoord &res, const GlobalCoordinate &coord, int dir)
Helper function for extracting given 2D coordinates from a 3D coordinate.
const Quad & getQuad(int col, int index) const
Definition: boundarygrid.hh:146
static BoundaryGrid uniform(const FaceCoord &min, const FaceCoord &max, int k1, int k2, bool dc=false)
Establish an uniform quad grid.
void add(const Quad &quad)
Add a quad to the grid.
bool find(Vertex &res, const Vertex &coord) const
Locate the position of a vertex on the grid.
Dune::FieldVector< double, 2 > FaceCoord
A coordinate on the quad grid.
Definition: boundarygrid.hh:38
size_t colSize(int i) const
Definition: boundarygrid.hh:157
size_t size() const
Obtain the number of quads in the grid.
Definition: boundarygrid.hh:152
Geometry class for general hexagons.
Definition: boundarygrid.hh:289
Isotropic linear elastic material.
Definition: materials.hh:26
A class for representing a general multi-point constraint equation.
Definition: mpc.hh:59
size_t getNoMaster() const
Returns the number of master DOFs.
Definition: mpc.hh:141
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:110
const DOF & getMaster(size_t i) const
Returns a reference to the i'th master DOF.
Definition: mpc.hh:138
static Material * create(int ID, const Dune::DynamicVector< double > ¶ms)
Creates a material object of a given type.
static Matrix augment(const Matrix &A, const Matrix &B, size_t r0, size_t c0, bool symmetric)
Augment a matrix with another.
static Matrix diagonal(size_t N)
Returns a diagonal matrix.
static Matrix fromDense(const Dune::DynamicMatrix< double > &T)
Create a sparse matrix from a dense matrix.
static void fromAdjacency(Matrix &A, const AdjacencyPattern &adj, int rows, int cols)
Create a sparse matrix from a given adjacency pattern.
Definition: mortar_schur.hpp:27
void apply(const Vector &x, Vector &y) const override
Apply the multiplier block.
Definition: mortar_schur.hpp:42
Definition: mortar_evaluator.hpp:27
Definition: mortar_schur_precond.hpp:37
Singleton handler for the set of LinearShapeFunction.
Definition: shapefunctions.hpp:181
@ n
Definition: shapefunctions.hpp:184
static const P1ShapeFunctionSet & instance()
Get the only instance of this class.
Definition: shapefunctions.hpp:193
Definition: shapefunctions.hpp:273
int size()
Definition: shapefunctions.hpp:322
Definition: uzawa_solver.hpp:25
Dune::MatrixAdapter< Matrix, Vector, Vector > Operator
A linear operator.
Definition: elasticity_preconditioners.hpp:49
Direction
An enum for specification of global coordinate directions.
Definition: mpc.hh:24
@ XYZ
Definition: mpc.hh:26
std::map< Material *, double > volume
Definition: elasticity_upscale_impl.hpp:659
upscaledRho
Definition: elasticity_upscale_impl.hpp:582
@ ITERATIVE
Definition: elasticity_upscale.hpp:54
BoundaryGrid::Vertex maxXmaxY(std::vector< BoundaryGrid::Vertex > &in)
Find the vertex in the vector with maximum X and maximum Y.
min[0]
Definition: elasticity_upscale_impl.hpp:146
Dune::BCRSMatrix< Dune::FieldMatrix< double, 1, 1 > > Matrix
A sparse matrix holding our operator.
Definition: matrixops.hpp:27
std::vector< int > cells
Definition: elasticity_upscale_impl.hpp:657
IMPL_FUNC(std::vector< BoundaryGrid::Vertex >, extractFace(Direction dir, ctype coord))
Definition: elasticity_upscale_impl.hpp:30
BoundaryGrid::Vertex minXminY(std::vector< BoundaryGrid::Vertex > &in)
Find the vertex in the vector with minimum X and minimum Y.
std::vector< double > Emod
Definition: elasticity_upscale_impl.hpp:578
int j
Definition: elasticity_upscale_impl.hpp:658
Dune::BlockVector< Dune::FieldVector< double, 1 > > Vector
A vector holding our RHS.
Definition: matrixops.hpp:33
const LeafVertexIterator itend
Definition: elasticity_upscale_impl.hpp:147
BoundaryGrid::Vertex maxXminY(std::vector< BoundaryGrid::Vertex > &in)
Find the vertex in the vector with maximum X and minimum Y.
static const int bfunc
Definition: elasticity_upscale_impl.hpp:444
int numsolvers
Definition: elasticity_upscale_impl.hpp:910
@ SCHURAMG
schur + amg
Definition: elasticity_upscale.hpp:69
@ SCHUR
schur + primary preconditioner
Definition: elasticity_upscale.hpp:68
@ SCHURTWOLEVEL
schur + twolevel
Definition: elasticity_upscale.hpp:71
@ SIMPLE
diagonal approximation of A
Definition: elasticity_upscale.hpp:67
@ SCHURSCHWARZ
schur + schwarz+lu
Definition: elasticity_upscale.hpp:70
master push_back(extractMasterFace(Y, min[1]))
std::vector< double > rho
Definition: elasticity_upscale_impl.hpp:581
MaterialMap cache
Definition: elasticity_upscale_impl.hpp:577
determineSideFaces(min, max)
std::vector< double > Poiss
Definition: elasticity_upscale_impl.hpp:579
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:193
Dune::FieldVector< ctype, comp > eps0
Definition: elasticity_upscale_impl.hpp:446
std::vector< int > satnum
Definition: elasticity_upscale_impl.hpp:580
BoundaryGrid::Vertex minXmaxY(std::vector< BoundaryGrid::Vertex > &in)
Find the vertex in the vector with minimum X and maximum Y.
@ AMG
Definition: elasticity_upscale.hpp:58
@ TWOLEVEL
Definition: elasticity_upscale.hpp:61
@ SCHWARZ
Definition: elasticity_upscale.hpp:60
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:892
auto deck
Definition: elasticity_upscale_impl.hpp:590
LeafVertexIterator start
Definition: elasticity_upscale_impl.hpp:150
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
Dune::Amg::AMG< Operator, Vector, Smoother > type
Definition: elasticity_preconditioners.hpp:106
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
Dune::Amg::TwoLevelMethod< Operator, CoarsePolicy, Schwarz::type > type
Definition: elasticity_preconditioners.hpp:169
int node
Node number identifying this DOF.
Definition: mpc.hh:80
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