OilfieldToolsNet/opm-upscaling
4
0
Fork 0
mirror of https://github.com/OPM/opm-upscaling.git synced 2025-02-25 18:45:23 -06:00
Code Issues Packages Projects Releases Wiki Activity

quell compiler warnings (sign/unsigned mismatch, unused variables)

Browse Source
This commit is contained in:
Arne Morten Kvarving 2012-11-20 11:22:04 +01:00
parent 6f1e511531
commit 343411cc9b
8 changed files with 42 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
dune/elasticity/asmhandler.hpp
View File

@ -55,7 +55,7 @@ class ASMHandler {
//! \brief Get the number of equations in the system //! \brief Get the number of equations in the system
//! \returns The number of equations in the system //! \returns The number of equations in the system
int getEqns() const size_t getEqns() const
{ {
return maxeqn; return maxeqn;
} }

20
dune/elasticity/asmhandler_impl.hpp
View File

@ -59,7 +59,7 @@ void ASMHandler<GridType>::addDOF(int row, int erow,
for (int l=0;l<dim;++l) { for (int l=0;l<dim;++l) {
MPC* mpc = getMPC(index2,l); MPC* mpc = getMPC(index2,l);
if (mpc) { if (mpc) {
for (int n=0;n<mpc->getNoMaster();++n) { for (size_t n=0;n<mpc->getNoMaster();++n) {
int idx = meqn[mpc->getMaster(n).node*dim+mpc->getMaster(n).dof-1]; int idx = meqn[mpc->getMaster(n).node*dim+mpc->getMaster(n).dof-1];
if (idx != -1) if (idx != -1)
A[row][idx] += scale*mpc->getMaster(n).coeff*(*K)[erow][j*dim+l]; A[row][idx] += scale*mpc->getMaster(n).coeff*(*K)[erow][j*dim+l];
@ -93,7 +93,7 @@ void ASMHandler<GridType>::addElement(
for (int k=0;k<dim;++k) { for (int k=0;k<dim;++k) {
MPC* mpc = getMPC(index1,k); MPC* mpc = getMPC(index1,k);
if (mpc) { if (mpc) {
for (int n=0;n<mpc->getNoMaster();++n) { for (size_t n=0;n<mpc->getNoMaster();++n) {
int idx = meqn[mpc->getMaster(n).node*dim+mpc->getMaster(n).dof-1]; int idx = meqn[mpc->getMaster(n).node*dim+mpc->getMaster(n).dof-1];
addDOF(idx,i*dim+k,K,S,set,cell,b2,mpc->getMaster(n).coeff); addDOF(idx,i*dim+k,K,S,set,cell,b2,mpc->getMaster(n).coeff);
} }
@ -124,7 +124,7 @@ void ASMHandler<GridType>::extractValues(Dune::FieldVector<double,comp>& v,
if (it2 != fixedNodes.end() && it2->second.first & (1 << n)) if (it2 != fixedNodes.end() && it2->second.first & (1 << n))
v[l++] = it2->second.second[n]; v[l++] = it2->second.second[n];
else if (mpc) { else if (mpc) {
for (int m=0;m<mpc->getNoMaster();++m) { for (size_t m=0;m<mpc->getNoMaster();++m) {
int idx = meqn[mpc->getMaster(m).node*dim+mpc->getMaster(m).dof-1]; int idx = meqn[mpc->getMaster(m).node*dim+mpc->getMaster(m).dof-1];
if (idx != -1) if (idx != -1)
v[l] += u[idx]*mpc->getMaster(m).coeff; v[l] += u[idx]*mpc->getMaster(m).coeff;
@ -143,7 +143,7 @@ void ASMHandler<GridType>::expandSolution(Vector& result, const Vector& u)
result.resize(nodes*dim); result.resize(nodes*dim);
result = 0; result = 0;
int l=0; int l=0;
for (size_t i=0;i<nodes;++i) { for (int i=0;i<nodes;++i) {
fixIt it = fixedNodes.find(i); fixIt it = fixedNodes.find(i);
Direction dir; Direction dir;
if (it == fixedNodes.end()) if (it == fixedNodes.end())
@ -163,11 +163,11 @@ void ASMHandler<GridType>::expandSolution(Vector& result, const Vector& u)
} }
// second loop - handle MPC couplings // second loop - handle MPC couplings
l = 0; l = 0;
for (size_t i=0;i<nodes;++i) { for (int i=0;i<nodes;++i) {
for (int j=0;j<dim;++j) { for (int j=0;j<dim;++j) {
MPC* mpc = getMPC(i,j); MPC* mpc = getMPC(i,j);
if (mpc) { if (mpc) {
for (int n=0;n<mpc->getNoMaster();++n) { for (size_t n=0;n<mpc->getNoMaster();++n) {
int idx = mpc->getMaster(n).node*dim+mpc->getMaster(n).dof-1; int idx = mpc->getMaster(n).node*dim+mpc->getMaster(n).dof-1;
if (meqn[idx] != -1) if (meqn[idx] != -1)
result[l] += u[meqn[idx]]*mpc->getMaster(n).coeff; result[l] += u[meqn[idx]]*mpc->getMaster(n).coeff;
@ -188,7 +188,7 @@ void ASMHandler<GridType>::addMPC(MPC* mpc)
fixIt it = fixedNodes.find(mpc->getSlave().node); fixIt it = fixedNodes.find(mpc->getSlave().node);
int flag = 1 << (mpc->getSlave().dof-1); int flag = 1 << (mpc->getSlave().dof-1);
if (it == fixedNodes.end() || if (it == fixedNodes.end() ||
!(it->second.first & (1 << (mpc->getSlave().dof-1)))) { !(it->second.first & flag)) {
mpcs.insert(std::make_pair(slaveNode,mpc)); mpcs.insert(std::make_pair(slaveNode,mpc));
return; return;
} }
@ -257,7 +257,7 @@ void ASMHandler<GridType>::resolveMPCChain(MPC* mpc)
coeff[master.dof-1] = master.coeff; coeff[master.dof-1] = master.coeff;
int removeOld = 0; int removeOld = 0;
for (size_t d = 1; d <= dim; d++) for (int d = 1; d <= dim; d++)
if (fabs(coeff[d-1]) > 1.0e-8) if (fabs(coeff[d-1]) > 1.0e-8)
{ {
MPC* mpc2 = getMPC(mpc->getMaster(i).node,d-1); MPC* mpc2 = getMPC(mpc->getMaster(i).node,d-1);
@ -354,7 +354,7 @@ void ASMHandler<GridType>::nodeAdjacency(const LeafIterator& it,
for (int l=0;l<dim;++l) { for (int l=0;l<dim;++l) {
MPC* mpc = getMPC(indexj,l); MPC* mpc = getMPC(indexj,l);
if (mpc) { if (mpc) {
for (int i=0;i<mpc->getNoMaster();++i) { for (size_t i=0;i<mpc->getNoMaster();++i) {
int idx = meqn[mpc->getMaster(i).node*dim+ int idx = meqn[mpc->getMaster(i).node*dim+
mpc->getMaster(i).dof-1]; mpc->getMaster(i).dof-1];
if (idx != -1) if (idx != -1)
@ -386,7 +386,7 @@ void ASMHandler<GridType>::determineAdjacencyPattern()
for (int k=0;k<dim;++k) { for (int k=0;k<dim;++k) {
MPC* mpc = getMPC(indexi,k); MPC* mpc = getMPC(indexi,k);
if (mpc) { if (mpc) {
for (int l=0;l<mpc->getNoMaster();++l) { for (size_t l=0;l<mpc->getNoMaster();++l) {
nodeAdjacency(it,vertexsize, nodeAdjacency(it,vertexsize,
meqn[mpc->getMaster(l).node*dim+ meqn[mpc->getMaster(l).node*dim+
mpc->getMaster(l).dof-1]); mpc->getMaster(l).dof-1]);

8
dune/elasticity/boundarygrid.cpp
View File

@ -175,7 +175,7 @@ int BoundaryGrid::Q4inv(FaceCoord& res, const Quad& q,
// check that obtained solutions are inside element // check that obtained solutions are inside element
double tol = 1+epsOut; double tol = 1+epsOut;
int nInside=0; int nInside=0;
for (int i=0;i<xi.size();++i) { for (size_t i=0;i<xi.size();++i) {
if (xi[i] < tol && eta[i] < tol) { if (xi[i] < tol && eta[i] < tol) {
if (++nInside > 1) { if (++nInside > 1) {
std::cout << "multiple solutions" << std::endl; std::cout << "multiple solutions" << std::endl;
@ -232,7 +232,7 @@ bool BoundaryGrid::bilinearSolve(double epsilon, double order,
double Q2 = A[0]*B[1]-B[0]*A[1]; double Q2 = A[0]*B[1]-B[0]*A[1];
std::vector<double> Z; std::vector<double> Z;
cubicSolve(epsilon,0,Q2,Q1,Q0,Z); cubicSolve(epsilon,0,Q2,Q1,Q0,Z);
for (int i=0;i<Z.size();++i) { for (size_t i=0;i<Z.size();++i) {
Q0 = A[0]*Z[i]+A[2]; Q0 = A[0]*Z[i]+A[2];
if (fabs(Q0) > tol) { if (fabs(Q0) > tol) {
X.push_back(Z[i]); X.push_back(Z[i]);
@ -244,10 +244,10 @@ bool BoundaryGrid::bilinearSolve(double epsilon, double order,
Q2 = A[0]*B[2]-B[0]*A[2]; Q2 = A[0]*B[2]-B[0]*A[2];
Z.clear(); Z.clear();
cubicSolve(epsilon,0,Q2,Q1,Q0,Z); cubicSolve(epsilon,0,Q2,Q1,Q0,Z);
for (int i=0;i<Z.size();++i) { for (size_t i=0;i<Z.size();++i) {
Q0 = A[0]*Z[i]+A[1]; Q0 = A[0]*Z[i]+A[1];
if (fabs(Q0) > tol) { if (fabs(Q0) > tol) {
int j=0; size_t j=0;
for (j=0;j<Y.size();++j) for (j=0;j<Y.size();++j)
if (fabs(Y[j]-Z[i]) <= epsilon*order) break; if (fabs(Y[j]-Z[i]) <= epsilon*order) break;
if (j == Y.size()) { if (j == Y.size()) {

5
dune/elasticity/boundarygrid.hh
View File

@ -237,7 +237,7 @@ class BoundaryGrid {
//! \brief Print to a stream //! \brief Print to a stream
friend std::ostream& operator <<(std::ostream& os, const BoundaryGrid& g) friend std::ostream& operator <<(std::ostream& os, const BoundaryGrid& g)
{ {
for (int i=0;i<g.size();++i) for (size_t i=0;i<g.size();++i)
os << g[i] << std::endl; os << g[i] << std::endl;
return os; return os;
} }
@ -326,7 +326,7 @@ class HexGeometry<2, cdim, GridImp>
const typename GridImp::LeafGridView::template Codim<3>::Iterator itend = gv.leafView().template end<3>(); const typename GridImp::LeafGridView::template Codim<3>::Iterator itend = gv.leafView().template end<3>();
for (int i=0;i<4;++i) { for (int i=0;i<4;++i) {
typename GridImp::LeafGridView::template Codim<3>::Iterator it=start; typename GridImp::LeafGridView::template Codim<3>::Iterator it=start;
for (it ; it != itend; ++it) { for (; it != itend; ++it) {
if (mapper.map(*it) == q.v[i].i) if (mapper.map(*it) == q.v[i].i)
break; break;
} }
@ -414,7 +414,6 @@ class HexGeometry<2, cdim, GridImp>
Dune::GenericReferenceElements< ctype, 2 >::general(type()); Dune::GenericReferenceElements< ctype, 2 >::general(type());
LocalCoordinate x = refElement.position(0,0); LocalCoordinate x = refElement.position(0,0);
LocalCoordinate dx; LocalCoordinate dx;
int i=0;
do { do {
using namespace Dune::GenericGeometry; using namespace Dune::GenericGeometry;
// DF^n dx^n = F^n, x^{n+1} -= dx^n // DF^n dx^n = F^n, x^{n+1} -= dx^n

6
dune/elasticity/elasticity_impl.hpp
View File

@ -30,7 +30,7 @@ void Elasticity<GridType>::getBmatrix(Dune::FieldMatrix<ctype,components,funcdim
if (dim == 3) { if (dim == 3) {
#define INDEX(i,j) i+6*j #define INDEX(i,j) i+6*j
for (size_t i=0; i < funcs; ++i) { for (int i=0; i < funcs; ++i) {
// normal strain part // normal strain part
temp[INDEX(0,0)][i] = dNdX[i][0]; temp[INDEX(0,0)][i] = dNdX[i][0];
temp[INDEX(1,1)][i] = dNdX[i][1]; temp[INDEX(1,1)][i] = dNdX[i][1];
@ -48,7 +48,7 @@ void Elasticity<GridType>::getBmatrix(Dune::FieldMatrix<ctype,components,funcdim
} else if (dim == 2) { } else if (dim == 2) {
#undef INDEX #undef INDEX
#define INDEX(i,j) i+3*j #define INDEX(i,j) i+3*j
for (size_t i=0; i < funcs; ++i) { for (int i=0; i < funcs; ++i) {
// normal strain part // normal strain part
temp[INDEX(0,0)][i] = dNdX[i][0]; temp[INDEX(0,0)][i] = dNdX[i][0];
temp[INDEX(1,1)][i] = dNdX[i][1]; temp[INDEX(1,1)][i] = dNdX[i][1];
@ -59,7 +59,7 @@ void Elasticity<GridType>::getBmatrix(Dune::FieldMatrix<ctype,components,funcdim
} }
} }
size_t k=0; size_t k=0;
for (size_t j=0;j<funcs*dim;++j) for (int j=0;j<funcs*dim;++j)
for (size_t i=0;i<components;++i,++k) for (size_t i=0;i<components;++i,++k)
B[i][j] = temp[k%rows][k/rows]; B[i][j] = temp[k%rows][k/rows];
} }

49
dune/elasticity/elasticity_upscale_impl.hpp
View File

@ -14,7 +14,6 @@
std::vector<BoundaryGrid::Vertex> ElasticityUpscale<GridType>::extractFace(Direction dir, ctype coord) std::vector<BoundaryGrid::Vertex> ElasticityUpscale<GridType>::extractFace(Direction dir, ctype coord)
{ {
std::vector<BoundaryGrid::Vertex> result; std::vector<BoundaryGrid::Vertex> result;
const LeafIndexSet& set = gv.leafView().indexSet();
const LeafVertexIterator itend = gv.leafView().template end<dim>(); const LeafVertexIterator itend = gv.leafView().template end<dim>();
// make a mapper for codim dim entities in the leaf grid // make a mapper for codim dim entities in the leaf grid
@ -68,7 +67,7 @@ BoundaryGrid ElasticityUpscale<GridType>::extractMasterFace(Direction dir,
else if (side == RIGHT) else if (side == RIGHT)
idx = set.subIndex(*cell,V2[i][0],dim); idx = set.subIndex(*cell,V2[i][0],dim);
LeafVertexIterator it=start; LeafVertexIterator it=start;
for (it ; it != itend; ++it) { for (; it != itend; ++it) {
if (mapper.map(*it) == idx) if (mapper.map(*it) == idx)
break; break;
} }
@ -79,7 +78,7 @@ BoundaryGrid ElasticityUpscale<GridType>::extractMasterFace(Direction dir,
if (side == RIGHT) if (side == RIGHT)
idx = set.subIndex(*cell,V2[i][j],dim); idx = set.subIndex(*cell,V2[i][j],dim);
LeafVertexIterator it=start; LeafVertexIterator it=start;
for (it ; it != itend; ++it) { for (; it != itend; ++it) {
if (mapper.map(*it) == idx) if (mapper.map(*it) == idx)
break; break;
} }
@ -145,7 +144,6 @@ void ElasticityUpscale<GridType>::findBoundaries(double* min,
{ {
max[0] = max[1] = max[2] = -1e5; max[0] = max[1] = max[2] = -1e5;
min[0] = min[1] = min[2] = 1e5; min[0] = min[1] = min[2] = 1e5;
const LeafIndexSet& set = gv.leafView().indexSet();
const LeafVertexIterator itend = gv.leafView().template end<dim>(); const LeafVertexIterator itend = gv.leafView().template end<dim>();
// iterate over vertices and find slaves // iterate over vertices and find slaves
@ -179,7 +177,7 @@ void ElasticityUpscale<GridType>::periodicPlane(Direction plane, Direction dir,
const std::vector<BoundaryGrid::Vertex>& slave, const std::vector<BoundaryGrid::Vertex>& slave,
const BoundaryGrid& master) const BoundaryGrid& master)
{ {
for (int i=0;i<slave.size();++i) { for (size_t i=0;i<slave.size();++i) {
BoundaryGrid::Vertex coord; BoundaryGrid::Vertex coord;
if (master.find(coord,slave[i])) { if (master.find(coord,slave[i])) {
addMPC(X,slave[i].i,coord); addMPC(X,slave[i].i,coord);
@ -204,7 +202,7 @@ Matrix ElasticityUpscale<GridType>::findBMatrixLLM(const SlaveGrid& slave)
continue; continue;
MPC* mpc = A.getMPC(it->first,k); MPC* mpc = A.getMPC(it->first,k);
if (mpc) { if (mpc) {
for (int n=0;n<mpc->getNoMaster();++n) { for (size_t n=0;n<mpc->getNoMaster();++n) {
int idx = A.getEquationForDof(mpc->getMaster(n).node, int idx = A.getEquationForDof(mpc->getMaster(n).node,
mpc->getMaster(n).dof-1); mpc->getMaster(n).dof-1);
if (idx != -1) if (idx != -1)
@ -227,7 +225,7 @@ Matrix ElasticityUpscale<GridType>::findBMatrixLLM(const SlaveGrid& slave)
continue; continue;
MPC* mpc = A.getMPC(it->first,k); MPC* mpc = A.getMPC(it->first,k);
if (mpc) { if (mpc) {
for (int n=0;n<mpc->getNoMaster();++n) { for (size_t n=0;n<mpc->getNoMaster();++n) {
int idx = A.getEquationForDof(mpc->getMaster(n).node, int idx = A.getEquationForDof(mpc->getMaster(n).node,
mpc->getMaster(n).dof-1); mpc->getMaster(n).dof-1);
if (idx != -1) if (idx != -1)
@ -251,7 +249,7 @@ Matrix ElasticityUpscale<GridType>::findLMatrixLLM(const SlaveGrid& slave,
int nbeqn=0; int nbeqn=0;
std::vector<int> dofmap(master.totalNodes()*dim,-1); std::vector<int> dofmap(master.totalNodes()*dim,-1);
int col=0; int col=0;
for (int i=0;i<master.totalNodes();++i) { for (size_t i=0;i<master.totalNodes();++i) {
if (master.isFixed(i)) { if (master.isFixed(i)) {
dofmap[i*dim ] = -1; dofmap[i*dim ] = -1;
dofmap[i*dim+1] = -1; dofmap[i*dim+1] = -1;
@ -541,9 +539,7 @@ void ElasticityUpscale<GridType>::assemble(int loadcase, bool matrix)
{ {
const int comp = 3+(dim-2)*3; const int comp = 3+(dim-2)*3;
static const int bfunc = 4+(dim-2)*4; static const int bfunc = 4+(dim-2)*4;
const int N = gv.size(dim);
const LeafIndexSet& set = gv.leafView().indexSet();
const LeafIterator itend = gv.leafView().template end<0>(); const LeafIterator itend = gv.leafView().template end<0>();
Dune::FieldMatrix<ctype,comp,comp> C; Dune::FieldMatrix<ctype,comp,comp> C;
@ -567,8 +563,6 @@ void ElasticityUpscale<GridType>::assemble(int loadcase, bool matrix)
materials[m++]->getConstitutiveMatrix(C); materials[m++]->getConstitutiveMatrix(C);
// determine geometry type of the current element and get the matching reference element // determine geometry type of the current element and get the matching reference element
Dune::GeometryType gt = it->type(); Dune::GeometryType gt = it->type();
const Dune::template GenericReferenceElement<ctype,dim> &ref =
Dune::GenericReferenceElements<ctype,dim>::general(gt);
Dune::FieldMatrix<ctype,dim*bfunc,dim*bfunc> Aq; Dune::FieldMatrix<ctype,dim*bfunc,dim*bfunc> Aq;
K = 0; K = 0;
@ -616,14 +610,9 @@ void ElasticityUpscale<GridType>::averageStress(Dune::FieldVector<ctype,comp>& s
return; return;
static const int bfunc = 4+(dim-2)*4; static const int bfunc = 4+(dim-2)*4;
const int N = gv.size(dim);
const LeafIndexSet& set = gv.leafView().indexSet();
const LeafIterator itend = gv.leafView().template end<0>(); const LeafIterator itend = gv.leafView().template end<0>();
// get a set of P1 shape functions
P1ShapeFunctionSet<ctype,ctype,dim> basis = P1ShapeFunctionSet<ctype,ctype,dim>::instance();
Dune::FieldMatrix<ctype,comp,comp> C; Dune::FieldMatrix<ctype,comp,comp> C;
Dune::FieldVector<ctype,comp> eps0; Dune::FieldVector<ctype,comp> eps0;
eps0 = 0; eps0 = 0;
@ -635,8 +624,6 @@ void ElasticityUpscale<GridType>::averageStress(Dune::FieldVector<ctype,comp>& s
materials[m++]->getConstitutiveMatrix(C); materials[m++]->getConstitutiveMatrix(C);
// determine geometry type of the current element and get the matching reference element // determine geometry type of the current element and get the matching reference element
Dune::GeometryType gt = it->type(); Dune::GeometryType gt = it->type();
const Dune::template GenericReferenceElement<ctype,dim> &ref =
Dune::GenericReferenceElements<ctype,dim>::general(gt);
Dune::FieldVector<ctype,bfunc*dim> v; Dune::FieldVector<ctype,bfunc*dim> v;
A.extractValues(v,u,it); A.extractValues(v,u,it);
@ -710,7 +697,7 @@ void ElasticityUpscale<GridType>::loadMaterialsFromGrid(const std::string& file)
// scale E modulus of materials // scale E modulus of materials
if (Escale > 0) { if (Escale > 0) {
Emin = *std::min_element(Emod.begin(),Emod.end()); Emin = *std::min_element(Emod.begin(),Emod.end());
for (int i=0;i<Emod.size();++i) for (size_t i=0;i<Emod.size();++i)
Emod[i] *= Escale/Emin; Emod[i] *= Escale/Emin;
} }
@ -770,9 +757,9 @@ void ElasticityUpscale<GridType>::loadMaterialsFromRocklist(const std::string& f
// scale E modulus of materials // scale E modulus of materials
if (Escale > 0) { if (Escale > 0) {
Emin=1e10; Emin=1e10;
for (int i=0;i<cache.size();++i) for (size_t i=0;i<cache.size();++i)
Emin = std::min(Emin,((Isotropic*)cache[i])->getE()); Emin = std::min(Emin,((Isotropic*)cache[i])->getE());
for (int i=0;i<cache.size();++i) { for (size_t i=0;i<cache.size();++i) {
double E = ((Isotropic*)cache[i])->getE(); double E = ((Isotropic*)cache[i])->getE();
((Isotropic*)cache[i])->setE(E*Escale/Emin); ((Isotropic*)cache[i])->setE(E*Escale/Emin);
} }
@ -780,7 +767,7 @@ void ElasticityUpscale<GridType>::loadMaterialsFromRocklist(const std::string& f
std::vector<double> volume; std::vector<double> volume;
volume.resize(cache.size()); volume.resize(cache.size());
if (file == "uniform") { if (file == "uniform") {
for (size_t i=0;i<gv.size(0);++i) for (int i=0;i<gv.size(0);++i)
materials.push_back(cache[0]); materials.push_back(cache[0]);
volume[0] = 1; volume[0] = 1;
} else { } else {
@ -951,8 +938,6 @@ void ElasticityUpscale<GridType>::periodicBCsLLM(const double* min,
Dune::LeafMultipleCodimMultipleGeomTypeMapper<GridType, Dune::LeafMultipleCodimMultipleGeomTypeMapper<GridType,
Dune::MCMGVertexLayout> mapper(gv); Dune::MCMGVertexLayout> mapper(gv);
LeafVertexIterator start=gv.leafView().template begin<dim>();
LeafVertexIterator itend = gv.leafView().template end<dim>();
BoundaryGrid::FaceCoord fmin,fmax; BoundaryGrid::FaceCoord fmin,fmax;
// YZ plane // YZ plane
@ -969,10 +954,10 @@ void ElasticityUpscale<GridType>::periodicBCsLLM(const double* min,
// step 5 // step 5
std::map<int,BoundaryGrid::Vertex> m; std::map<int,BoundaryGrid::Vertex> m;
// find matching coefficients // find matching coefficients
for (int i=0;i<master.size();++i) { for (size_t i=0;i<master.size();++i) {
for (int j=0;j<2;++j) { for (int j=0;j<2;++j) {
m.clear(); m.clear();
for (int k=0;k<slave[i*2+j].size();++k) { for (size_t k=0;k<slave[i*2+j].size();++k) {
BoundaryGrid::Vertex c; BoundaryGrid::Vertex c;
if (master[i].find(c,slave[i*2+j][k])) { if (master[i].find(c,slave[i*2+j][k])) {
m.insert(std::make_pair(slave[i*2+j][k].i,c)); m.insert(std::make_pair(slave[i*2+j][k].i,c));
@ -999,12 +984,12 @@ void ElasticityUpscale<GridType>::setupPreconditioner()
if (cell / gv.logicalCartesianSize()[2] > 0 if (cell / gv.logicalCartesianSize()[2] > 0
&& cell % gv.logicalCartesianSize()[2] == 0) && cell % gv.logicalCartesianSize()[2] == 0)
currdomain++; currdomain++;
for (int i=0;i<8;++i) { for (size_t i=0;i<8;++i) {
int idx=set.subIndex(*it,i,dim); int idx=set.subIndex(*it,i,dim);
for (int d=0;d<3;++d) { for (int d=0;d<3;++d) {
MPC* mpc = A.getMPC(idx,d); MPC* mpc = A.getMPC(idx,d);
if (mpc) { if (mpc) {
for (int n=0;n<mpc->getNoMaster();++n) { for (size_t n=0;n<mpc->getNoMaster();++n) {
int row = A.getEquationForDof(mpc->getMaster(n).node,d); int row = A.getEquationForDof(mpc->getMaster(n).node,d);
if (row > -1) if (row > -1)
rows[currdomain].insert(row); rows[currdomain].insert(row);
@ -1034,11 +1019,11 @@ void ElasticityUpscale<GridType>::setupSolvers(Solver solver)
// [ 0 0 0 L3' L4' 0 0] [ub_2] [0] // [ 0 0 0 L3' L4' 0 0] [ub_2] [0]
int r = A.getOperator().N(); int r = A.getOperator().N();
int c = A.getOperator().M(); int c = A.getOperator().M();
for (int i=0;i<B.size();++i) { for (size_t i=0;i<B.size();++i) {
A.getOperator() = MatrixOps::augment(A.getOperator(),B[i],0,c,true); A.getOperator() = MatrixOps::augment(A.getOperator(),B[i],0,c,true);
c += B[i].M(); c += B[i].M();
} }
for (int i=0;i<L.size();++i) { for (size_t i=0;i<L.size();++i) {
A.getOperator() = MatrixOps::augment(A.getOperator(),L[i],r,c,true); A.getOperator() = MatrixOps::augment(A.getOperator(),L[i],r,c,true);
r += L[i].N(); r += L[i].N();
if (i % 2 == 1) if (i % 2 == 1)
@ -1055,7 +1040,7 @@ void ElasticityUpscale<GridType>::setupSolvers(Solver solver)
// [L1' 0 0] [l_1] = [0] // [L1' 0 0] [l_1] = [0]
// [L2' 0 0] [l_2] [0] // [L2' 0 0] [l_2] [0]
int c = A.getOperator().M(); int c = A.getOperator().M();
for (int i=0;i<L.size();++i) { for (size_t i=0;i<L.size();++i) {
A.getOperator() = MatrixOps::augment(A.getOperator(),L[i],0,c,true); A.getOperator() = MatrixOps::augment(A.getOperator(),L[i],0,c,true);
c += L[i].M(); c += L[i].M();
} }

6
dune/elasticity/matrixops_impl.hpp
View File

@ -19,7 +19,7 @@ void MatrixOps::fromAdjacency(Matrix& A, const std::vector< std::set<int> >& adj
A.setSize(rows, cols, sum); A.setSize(rows, cols, sum);
A.setBuildMode(Matrix::random); A.setBuildMode(Matrix::random);
for (int i = 0; i < adj.size(); i++) for (size_t i = 0; i < adj.size(); i++)
A.setrowsize(i,adj[i].size()); A.setrowsize(i,adj[i].size());
A.endrowsizes(); A.endrowsizes();
@ -53,7 +53,7 @@ Matrix MatrixOps::Axpy(const Matrix& A, const Matrix& B, double alpha)
assert(A.M() == B.M() && A.N() == B.N()); assert(A.M() == B.M() && A.N() == B.N());
// establish union adjacency pattern // establish union adjacency pattern
std::vector<std::set<int> > adj; AdjacencyPattern adj;
adj.resize(A.N()); adj.resize(A.N());
for (Matrix::ConstIterator it = A.begin(); for (Matrix::ConstIterator it = A.begin();
it != A.end(); ++it) { it != A.end(); ++it) {
@ -120,7 +120,7 @@ Matrix MatrixOps::augment(const Matrix& A, const Matrix& B,
} }
if (symmetric) { if (symmetric) {
// always establish diagonal elements or superLU crashes // always establish diagonal elements or superLU crashes
for (int i=0;i<nrow;++i) for (size_t i=0;i<nrow;++i)
adj[i].insert(i); adj[i].insert(i);
} }
Matrix result; Matrix result;

5
examples/upscale_elasticity.cpp
View File

@ -205,7 +205,7 @@ void writeOutput(const Params& p, Opm::time::StopWatch& watch, int cells,
} }
f << "#" << std::endl f << "#" << std::endl
<<"# Materials: " << volume.size() << std::endl; <<"# Materials: " << volume.size() << std::endl;
for (int i=0;i<volume.size();++i) for (size_t i=0;i<volume.size();++i)
f << "#\t Material" << i+1 << ": " << volume[i]*100 << "%" << std::endl; f << "#\t Material" << i+1 << ": " << volume[i]*100 << "%" << std::endl;
f << "#" << std::endl f << "#" << std::endl
<< "######################################################################" << std::endl << "######################################################################" << std::endl
@ -217,7 +217,6 @@ int main(int argc, char** argv)
{ {
try { try {
static const int dim = 3; static const int dim = 3;
static const int bfunc = 4+(dim-2)*4;
typedef Dune::CpGrid GridType; typedef Dune::CpGrid GridType;
@ -287,7 +286,7 @@ int main(int argc, char** argv)
upscale.solve(p.solver,p.ltol,i); upscale.solve(p.solver,p.ltol,i);
upscale.A.expandSolution(field[i],upscale.u[i]); upscale.A.expandSolution(field[i],upscale.u[i]);
#define CLAMP(x) (fabs(x)<1.e-5?0.f:x) #define CLAMP(x) (fabs(x)<1.e-5?0.f:x)
for (int j=0;j<field[i].size();++j) { for (size_t j=0;j<field[i].size();++j) {
double val = field[i][j]; double val = field[i][j];
field[i][j] = CLAMP(val); field[i][j] = CLAMP(val);
} }