Loading...
Searching...
No Matches
asmhandler_impl.hpp
Go to the documentation of this file.
//==============================================================================
//==============================================================================
#ifndef ASMHANDLER_IMPL_HPP_
#define ASMHANDLER_IMPL_HPP_
#include <dune/common/version.hh>
#include <iostream>
namespace Opm {
namespace Elasticity {
template<class GridType>
void ASMHandler<GridType>::initForAssembly()
{
resolveMPCChains();
preprocess();
determineAdjacencyPattern();
MatrixOps::fromAdjacency(A,adjacencyPattern,
adjacencyPattern.size(),adjacencyPattern.size());
b.resize(adjacencyPattern.size());
b = 0;
adjacencyPattern.clear();
// print some information
std::cout << "\tNumber of nodes: " << gv.size(dim) << std::endl;
std::cout << "\tNumber of elements: " << gv.size(0) << std::endl;
std::cout << "\tNumber of constraints: " << mpcs.size() << std::endl;
int fixedDofs=0;
for (fixIt it = fixedNodes.begin(); it != fixedNodes.end(); ++it) {
if (it->second.first & X)
fixedDofs++;
if (it->second.first & Y)
fixedDofs++;
if (it->second.first & Z)
fixedDofs++;
}
std::cout << "\tNumber of fixed dofs: " << fixedDofs << std::endl;
}
template<class GridType>
template<int esize>
void ASMHandler<GridType>::addDOF(int row, int erow,
const Dune::FieldMatrix<double,esize,esize>* K,
const Dune::FieldVector<double,esize>* S,
const LeafIndexSet& set,
const LeafIterator& cell,
Vector* bptr,
double scale)
{
if (row == -1)
return;
if (K) {
for (int j=0;j<esize/dim;++j) {
int index2 = set.subIndex(*cell,j,dim);
for (int l=0;l<dim;++l) {
MPC* mpc = getMPC(index2,l);
if (mpc) {
for (size_t n=0;n<mpc->getNoMaster();++n) {
int idx = meqn[mpc->getMaster(n).node*dim+mpc->getMaster(n).dof-1];
if (idx != -1)
A[row][idx] += scale*mpc->getMaster(n).coeff*(*K)[erow][j*dim+l];
}
} else if (meqn[index2*dim+l] != -1) {
A[row][meqn[index2*dim+l]] += scale*(*K)[erow][j*dim+l];
}
}
}
}
if (S && bptr)
(*bptr)[row] += scale*(*S)[erow];
}
template<class GridType>
template<int esize>
void ASMHandler<GridType>::addElement(
const Dune::FieldMatrix<double,esize,esize>* K,
const Dune::FieldVector<double,esize>* S,
const LeafIterator& cell,
Vector* b2)
{
if (!b2)
b2 = &b;
const LeafIndexSet& set = gv.leafGridView().indexSet();
for (int i=0;i<esize/dim;++i) {
int index1 = set.subIndex(*cell,i,dim);
fixIt it = fixedNodes.find(index1);
if (it != fixedNodes.end() && it->second.first == XYZ)
continue;
for (int k=0;k<dim;++k) {
MPC* mpc = getMPC(index1,k);
if (mpc) {
for (size_t n=0;n<mpc->getNoMaster();++n) {
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);
}
} else
addDOF(meqn[index1*dim+k],i*dim+k,K,S,set,cell,b2);
}
}
}
template<class GridType>
template<int comp>
void ASMHandler<GridType>::extractValues(Dune::FieldVector<double,comp>& v,
const Vector& u,
const LeafIterator& it)
{
v = 0;
const LeafIndexSet& set = gv.leafGridView().indexSet();
auto ref = Dune::ReferenceElements<double,dim>::cube();
int vertexsize = ref.size(dim);
int l=0;
for (int i=0;i<vertexsize;++i) {
int indexi = set.subIndex(*it,i,dim);
fixIt it2 = fixedNodes.find(indexi);
for (int n=0;n<dim;++n) {
MPC* mpc = getMPC(indexi,n);
if (it2 != fixedNodes.end() && it2->second.first & (1 << n))
v[l++] = it2->second.second[n];
else if (mpc) {
for (size_t m=0;m<mpc->getNoMaster();++m) {
int idx = meqn[mpc->getMaster(m).node*dim+mpc->getMaster(m).dof-1];
if (idx != -1)
v[l] += u[idx]*mpc->getMaster(m).coeff;
}
l++;
} else
v[l++] = u[meqn[indexi*dim+n]];
}
}
}
template<class GridType>
void ASMHandler<GridType>::expandSolution(Vector& result, const Vector& u)
{
int nodes = gv.size(dim);
result.resize(nodes*dim);
result = 0;
int l=0;
for (int i=0;i<nodes;++i) {
fixIt it = fixedNodes.find(i);
Direction dir;
if (it == fixedNodes.end())
dir = NONE;
else
dir = it->second.first;
int flag=1;
for (int j=0;j<dim;++j) {
if (dir & flag)
result[l] = it->second.second[j];
else if (meqn[l] != -1)
result[l] = u[meqn[l]];
l++;
flag *= 2;
}
}
// second loop - handle MPC couplings
l = 0;
for (int i=0;i<nodes;++i) {
for (int j=0;j<dim;++j) {
const MPC* mpc = getMPC(i,j);
if (mpc) {
for (size_t n=0;n<mpc->getNoMaster();++n) {
int idx = mpc->getMaster(n).node*dim+mpc->getMaster(n).dof-1;
if (meqn[idx] != -1)
result[l] += u[meqn[idx]]*mpc->getMaster(n).coeff;
}
}
++l;
}
}
}
template<class GridType>
void ASMHandler<GridType>::addMPC(MPC* mpc)
{
if (!mpc)
return;
int slaveNode = mpc->getSlave().node*dim+mpc->getSlave().dof-1;
fixIt it = fixedNodes.find(mpc->getSlave().node);
int flag = 1 << (mpc->getSlave().dof-1);
if (it == fixedNodes.end() ||
!(it->second.first & flag)) {
mpcs.insert(std::make_pair(slaveNode,mpc));
return;
}
delete mpc;
}
template<class GridType>
MPC* ASMHandler<GridType>::getMPC(int node, int dof)
{
if (mpcs.find(node*dim+dof) != mpcs.end())
return mpcs[node*dim+dof];
return NULL;
}
template<class GridType>
void ASMHandler<GridType>::updateFixedNode(int node,
const std::pair<Direction,NodeValue>& entry)
{
fixIt it = fixedNodes.find(node);
// same type or new - update values/add and return
if (it == fixedNodes.end() || it->second.first == entry.first) {
fixedNodes[node] = entry;
return;
}
int temp = it->second.first;
temp |= entry.first;
it->second.first = (Direction)temp;
int flag = 1;
for (int i=0;i<dim;++i) {
if (entry.first & flag)
it->second.second[i] = entry.second[i];
flag *= 2;
}
}
template<class GridType>
void ASMHandler<GridType>::printOperator() const
{
MatrixOps::print(A);
}
template<class GridType>
void ASMHandler<GridType>::printLoadVector() const
{
for (int i=0;i<b.size();++i) {
double val = b[i];
std::cout << (fabs(val)>1.e-12?val:0.f) << std::endl;
}
}
template<class GridType>
void ASMHandler<GridType>::resolveMPCChain(MPC* mpc)
{
size_t nMaster = mpc->getNoMaster();
if (nMaster == 0) return; // no masters, prescribed displacement only
for (size_t i = 0; i < nMaster; i++)
{
// Check if the master node has a local coordinate system attached. If yes,
// the slave DOF might be coupled to all (local) DOFs of the master node.
const MPC::DOF& master = mpc->getMaster(i);
Dune::FieldVector<double,dim> coeff;
coeff = 0;
coeff[master.dof-1] = master.coeff;
int removeOld = 0;
for (int d = 1; d <= dim; d++)
if (fabs(coeff[d-1]) > 1.0e-8)
{
MPC* mpc2 = getMPC(mpc->getMaster(i).node,d-1);
if (mpc2)
{
// We have a master DOF which is a slave in another MPC.
// Invoke resolveMPCchain recursively to ensure that all master DOFs
// of that equation are not slaves themselves.
resolveMPCChain(mpc2);
// Remove current master specification, unless it has been updated
if (!removeOld) removeOld = 1;
// Add constant offset from the other equation
mpc->addOffset(mpc2->getSlave().coeff*coeff[d-1]);
// Add masters from the other equations
for (size_t j = 0; j < mpc2->getNoMaster(); j++)
mpc->addMaster(mpc2->getMaster(j).node,
mpc2->getMaster(j).dof,
mpc2->getMaster(j).coeff*coeff[d-1]);
}
else
// The master node is free, but has a local coordinate system
if (d != mpc->getMaster(i).dof)
// Add coupling to the other local DOFs of this master node.
mpc->addMaster(mpc->getMaster(i).node,d,coeff[d-1]);
else if (coeff[d-1] != mpc->getMaster(i).coeff)
{
// Update the coupling coefficient of this master DOF
// due to the local-to-global transformation
mpc->updateMaster(i,coeff[d-1]);
removeOld = -1;
}
}
else if (d == mpc->getMaster(i).dof && !removeOld)
// The coefficient of the current master DOF is zero,
removeOld = 1; // so remove it from the contraint equation
if (removeOld == 1)
{
// Remove the old master DOF specification when it has been replaced
mpc->removeMaster(i--);
nMaster--; // we don't need to check the added masters
}
}
}
template<class GridType>
void ASMHandler<GridType>::preprocess()
{
int nodes = gv.size(dim);
meqn.resize(nodes*dim);
// iterate over nodes
for (int indexi=0;indexi<nodes;++indexi) {
fixIt it2 = fixedNodes.find(indexi);
if (it2 == fixedNodes.end()) {
for (int i=0;i<dim;++i) {
const MPC* mpc = getMPC(indexi,i);
if (!mpc)
meqn[indexi*dim+i] = maxeqn++;
else
meqn[indexi*dim+i] = -1;
}
} else {
int flag=1;
for (int i=0;i<dim;++i) {
if (it2->second.first & flag)
meqn[indexi*dim+i] = -1;
else {
const MPC* mpc = getMPC(indexi,i);
if (!mpc)
meqn[indexi*dim+i] = maxeqn++;
else
meqn[indexi*dim+i] = -1;
}
flag *= 2;
}
}
}
std::cout << "\tnumber of equations: " << maxeqn << std::endl;
}
template<class GridType>
void ASMHandler<GridType>::nodeAdjacency(const LeafIterator& it,
int vertexsize, int row)
{
if (row == -1)
return;
const LeafIndexSet& set = gv.leafGridView().indexSet();
for (int j=0;j<vertexsize;++j) {
int indexj = set.subIndex(*it,j,dim);
for (int l=0;l<dim;++l) {
MPC* mpc = getMPC(indexj,l);
if (mpc) {
for (size_t i=0;i<mpc->getNoMaster();++i) {
int idx = meqn[mpc->getMaster(i).node*dim+
mpc->getMaster(i).dof-1];
if (idx != -1)
adjacencyPattern[row].insert(idx);
}
} else if (meqn[indexj*dim+l] != -1)
adjacencyPattern[row].insert(meqn[indexj*dim+l]);
}
}
}
template<class GridType>
void ASMHandler<GridType>::determineAdjacencyPattern()
{
adjacencyPattern.resize(maxeqn);
std::cout << "\tsetting up sparsity pattern..." << std::endl;
LoggerHelper help(gv.size(0), 5, 50000);
const LeafIndexSet& set = gv.leafGridView().indexSet();
LeafIterator itend = gv.leafGridView().template end<0>();
// iterate over cells
int cell=0;
for (LeafIterator it = gv.leafGridView().template begin<0>(); it != itend; ++it, ++cell) {
auto ref = Dune::ReferenceElements<double,dim>::cube();
int vertexsize = ref.size(dim);
for (int i=0; i < vertexsize; i++) {
int indexi = set.subIndex(*it,i,dim);
for (int k=0;k<dim;++k) {
const MPC* mpc = getMPC(indexi,k);
if (mpc) {
for (size_t l=0;l<mpc->getNoMaster();++l) {
nodeAdjacency(it,vertexsize,
meqn[mpc->getMaster(l).node*dim+
mpc->getMaster(l).dof-1]);
}
} else
nodeAdjacency(it,vertexsize,meqn[indexi*dim+k]);
}
}
if (cell % 10000 == 0)
help.log(cell, "\t\t... still processing ... cell ");
}
}
}} // namespace Opm, Elasticity
#endif
LoggerHelper
Helper class for progress logging during time consuming processes.
Definition logutils.hpp:21
LoggerHelper::log
void log(int it, const std::string &prefix)
Log to the terminal.
Definition logutils.hpp:54
Opm::Elasticity::ASMHandler::dim
static const int dim
The dimension of the grid.
Definition asmhandler.hpp:38
Opm::Elasticity::ASMHandler::initForAssembly
void initForAssembly()
This function needs to be called before starting the element assembly.
Definition asmhandler_impl.hpp:22
Opm::Elasticity::ASMHandler::mpcs
MPCMap mpcs
The set of MPC.
Definition asmhandler.hpp:210
Opm::Elasticity::ASMHandler::resolveMPCChain
void resolveMPCChain(MPC *mpc)
Internal function.
Definition asmhandler_impl.hpp:250
Opm::Elasticity::ASMHandler::gv
const GridType & gv
A reference to the grid in use.
Definition asmhandler.hpp:237
Opm::Elasticity::ASMHandler::fixedNodes
fixMap fixedNodes
The map holding information about our fixed nodes.
Definition asmhandler.hpp:225
Opm::Elasticity::ASMHandler::maxeqn
size_t maxeqn
The number of equations in the system.
Definition asmhandler.hpp:240
Opm::Elasticity::ASMHandler::adjacencyPattern
AdjacencyPattern adjacencyPattern
Holds the adjacency pattern of the sparse matrix.
Definition asmhandler.hpp:228
Opm::Elasticity::ASMHandler::meqn
std::vector< int > meqn
Vector of (interleaved) dof<->eqn mapping.
Definition asmhandler.hpp:213
Opm::Elasticity::ASMHandler::fixIt
fixMap::iterator fixIt
Iterator over a fixmap.
Definition asmhandler.hpp:222
Opm::Elasticity::ASMHandler::LeafIndexSet
GridType::LeafGridView::IndexSet LeafIndexSet
A set of indices.
Definition asmhandler.hpp:41
Opm::Elasticity::ASMHandler::addElement
void addElement(const Dune::FieldMatrix< double, esize, esize > *K, const Dune::FieldVector< double, esize > *S, const LeafIterator &cell, Vector *b=NULL)
Add an element matrix/vector to the system.
Definition asmhandler_impl.hpp:85
Opm::Elasticity::ASMHandler::resolveMPCChains
void resolveMPCChains()
Resolve chained MPCs.
Definition asmhandler.hpp:168
Opm::Elasticity::ASMHandler::b
Vector b
The load vector.
Definition asmhandler.hpp:234
Opm::Elasticity::ASMHandler::A
Matrix A
The linear operator.
Definition asmhandler.hpp:231
Opm::Elasticity::ASMHandler::updateFixedNode
void updateFixedNode(int node, const std::pair< Direction, NodeValue > &entry)
Update/add a fixed node.
Definition asmhandler_impl.hpp:214
Opm::Elasticity::ASMHandler::extractValues
void extractValues(Dune::FieldVector< double, comp > &v, const Vector &u, const LeafIterator &it)
Extract values corresponding to cell.
Definition asmhandler_impl.hpp:114
Opm::Elasticity::ASMHandler::expandSolution
void expandSolution(Vector &result, const Vector &u)
Expand a system vector to a solution vector.
Definition asmhandler_impl.hpp:145
Opm::Elasticity::ASMHandler::addMPC
void addMPC(MPC *mpc)
Add a MPC.
Definition asmhandler_impl.hpp:187
Opm::Elasticity::ASMHandler::preprocess
void preprocess()
Internal function. Generate meqn for registered MPC/fixed nodes.
Definition asmhandler_impl.hpp:315
Opm::Elasticity::ASMHandler::nodeAdjacency
void nodeAdjacency(const LeafIterator &it, int vertexsize, int row)
Internal function.
Definition asmhandler_impl.hpp:351
Opm::Elasticity::ASMHandler::LeafIterator
GridType::LeafGridView::template Codim< 0 >::Iterator LeafIterator
An iterator over grid cells.
Definition asmhandler.hpp:44
Opm::Elasticity::ASMHandler::printLoadVector
void printLoadVector() const
Print the current load vector.
Definition asmhandler_impl.hpp:241
Opm::Elasticity::ASMHandler::getMPC
MPC * getMPC(int node, int dof)
Look for and return a MPC for a specified node+dof.
Definition asmhandler_impl.hpp:205
Opm::Elasticity::ASMHandler::determineAdjacencyPattern
void determineAdjacencyPattern()
Internal function. Calculate adjacency pattern.
Definition asmhandler_impl.hpp:375
Opm::Elasticity::ASMHandler::addDOF
void addDOF(int row, int erow, const Dune::FieldMatrix< double, esize, esize > *K, const Dune::FieldVector< double, esize > *S, const LeafIndexSet &set, const LeafIterator &cell, Vector *b, double scale=1.f)
Internal function.
Definition asmhandler_impl.hpp:52
Opm::Elasticity::ASMHandler::printOperator
void printOperator() const
Print the current operator.
Definition asmhandler_impl.hpp:235
Opm::Elasticity::MPC
A class for representing a general multi-point constraint equation.
Definition mpc.hh:59
Opm::Elasticity::MPC::getNoMaster
size_t getNoMaster() const
Returns the number of master DOFs.
Definition mpc.hh:141
Opm::Elasticity::MPC::addMaster
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
Opm::Elasticity::MPC::getSlave
const DOF & getSlave() const
Returns a reference to the slave DOF.
Definition mpc.hh:136
Opm::Elasticity::MPC::addOffset
void addOffset(double offset)
Increments the c0 coefficient by a given offset.
Definition mpc.hh:130
Opm::Elasticity::MPC::removeMaster
void removeMaster(size_t pos)
Removes the pos'th master DOF from the constraint equation.
Definition mpc.hh:123
Opm::Elasticity::MPC::updateMaster
void updateMaster(size_t pos, double c)
Updates the coefficient of the pos'th master DOF.
Definition mpc.hh:116
Opm::Elasticity::MPC::getMaster
const DOF & getMaster(size_t i) const
Returns a reference to the i'th master DOF.
Definition mpc.hh:138
Opm::Elasticity::MatrixOps::print
static void print(const Matrix &A)
Print a matrix to stdout.
Definition matrixops.cpp:73
Opm::Elasticity::MatrixOps::fromAdjacency
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
Opm::Elasticity::Vector
Dune::BlockVector< Dune::FieldVector< double, 1 > > Vector
A vector holding our RHS.
Definition matrixops.hpp:33
Opm::Elasticity::Direction
Direction
An enum for specification of global coordinate directions.
Definition mpc.hh:24
Opm
Inverting small matrices.
Definition ImplicitAssembly.hpp:43
Opm::Elasticity::MPC::DOF
A struct for representing one term (DOF number and associated coefficient) in a MPC equation.
Definition mpc.hh:67
Opm::Elasticity::MPC::DOF::dof
int dof
Local DOF number within node.
Definition mpc.hh:81
Opm::Elasticity::MPC::DOF::node
int node
Node number identifying this DOF.
Definition mpc.hh:80
Opm::Elasticity::MPC::DOF::coeff
double coeff
The constrained value, or master DOF scaling coefficient.
Definition mpc.hh:82