12 #ifndef ASMHANDLER_IMPL_HPP_
13 #define ASMHANDLER_IMPL_HPP_
15 #include <dune/common/version.hh>
19 namespace Elasticity {
21 template<
class Gr
idType>
26 determineAdjacencyPattern();
29 adjacencyPattern.size(),adjacencyPattern.size());
30 b.resize(adjacencyPattern.size());
32 adjacencyPattern.clear();
35 std::cout <<
"\tNumber of nodes: " << gv.size(dim) << std::endl;
36 std::cout <<
"\tNumber of elements: " << gv.size(0) << std::endl;
37 std::cout <<
"\tNumber of constraints: " << mpcs.size() << std::endl;
39 for (
fixIt it = fixedNodes.begin(); it != fixedNodes.end(); ++it) {
40 if (it->second.first &
X)
42 if (it->second.first &
Y)
44 if (it->second.first &
Z)
47 std::cout <<
"\tNumber of fixed dofs: " << fixedDofs << std::endl;
50 template<
class Gr
idType>
53 const Dune::FieldMatrix<double,esize,esize>* K,
54 const Dune::FieldVector<double,esize>* S,
63 for (
int j=0;
j<esize/dim;++
j) {
64 int index2 = set.subIndex(*cell,
j,dim);
65 for (
int l=0;l<dim;++l) {
66 MPC* mpc = getMPC(index2,l);
73 }
else if (meqn[index2*dim+l] != -1) {
74 A[row][meqn[index2*dim+l]] += scale*(*K)[erow][
j*dim+l];
80 (*b)[row] += scale*(*S)[erow];
83 template<
class Gr
idType>
86 const Dune::FieldMatrix<double,esize,esize>* K,
87 const Dune::FieldVector<double,esize>* S,
94 for (
int i=0;i<esize/dim;++i) {
95 int index1 = set.subIndex(*cell,i,dim);
96 fixIt it = fixedNodes.find(index1);
97 if (it != fixedNodes.end() && it->second.first ==
XYZ)
99 for (
int k=0;k<dim;++k) {
100 MPC* mpc = getMPC(index1,k);
107 addDOF(meqn[index1*dim+k],i*dim+k,K,S,set,cell,b2);
112 template<
class Gr
idType>
120 Dune::GeometryType gt = it->type();
122 const Dune::template ReferenceElement<double,dim> &ref =
123 Dune::ReferenceElements<double,dim>::general(gt);
124 int vertexsize = ref.size(dim);
126 for (
int i=0;i<vertexsize;++i) {
127 int indexi = set.subIndex(*it,i,dim);
128 fixIt it2 = fixedNodes.find(indexi);
129 for (
int n=0;n<dim;++n) {
130 MPC* mpc = getMPC(indexi,n);
131 if (it2 != fixedNodes.end() && it2->second.first & (1 << n))
132 v[l++] = it2->second.second[n];
141 v[l++] = u[meqn[indexi*dim+n]];
146 template<
class Gr
idType>
149 int nodes = gv.size(dim);
150 result.resize(nodes*dim);
153 for (
int i=0;i<nodes;++i) {
154 fixIt it = fixedNodes.find(i);
156 if (it == fixedNodes.end())
159 dir = it->second.first;
162 for (
int j=0;
j<dim;++
j) {
164 result[l] = it->second.second[
j];
165 else if (meqn[l] != -1)
166 result[l] = u[meqn[l]];
173 for (
int i=0;i<nodes;++i) {
174 for (
int j=0;
j<dim;++
j) {
175 MPC* mpc = getMPC(i,
j);
188 template<
class Gr
idType>
197 if (it == fixedNodes.end() ||
198 !(it->second.first & flag)) {
199 mpcs.insert(std::make_pair(slaveNode,mpc));
206 template<
class Gr
idType>
209 if (mpcs.find(node*dim+dof) != mpcs.end())
210 return mpcs[node*dim+dof];
215 template<
class Gr
idType>
217 const std::pair<Direction,NodeValue>& entry)
219 fixIt it = fixedNodes.find(node);
221 if (it == fixedNodes.end() || it->second.first == entry.first) {
222 fixedNodes[node] = entry;
225 int temp = it->second.first;
229 for (
int i=0;i<dim;++i) {
230 if (entry.first & flag)
231 it->second.second[i] = entry.second[i];
236 template<
class Gr
idType>
242 template<
class Gr
idType>
245 for (
int i=0;i<b.size();++i) {
247 std::cout << (fabs(val)>1.e-12?val:0.f) << std::endl;
251 template<
class Gr
idType>
255 if (nMaster == 0)
return;
257 for (
size_t i = 0; i < nMaster; i++)
262 Dune::FieldVector<double,dim> coeff;
267 for (
int d = 1; d <= dim; d++)
268 if (fabs(coeff[d-1]) > 1.0e-8)
276 resolveMPCChain(mpc2);
279 if (!removeOld) removeOld = 1;
316 template<
class Gr
idType>
319 int nodes = gv.size(dim);
320 meqn.resize(nodes*dim);
323 for (
int indexi=0;indexi<nodes;++indexi) {
324 fixIt it2 = fixedNodes.find(indexi);
325 if (it2 == fixedNodes.end()) {
326 for (
int i=0;i<dim;++i) {
327 MPC* mpc = getMPC(indexi,i);
329 meqn[indexi*dim+i] = maxeqn++;
331 meqn[indexi*dim+i] = -1;
335 for (
int i=0;i<dim;++i) {
336 if (it2->second.first & flag)
337 meqn[indexi*dim+i] = -1;
339 MPC* mpc = getMPC(indexi,i);
341 meqn[indexi*dim+i] = maxeqn++;
343 meqn[indexi*dim+i] = -1;
349 std::cout <<
"\tnumber of equations: " << maxeqn << std::endl;
352 template<
class Gr
idType>
354 int vertexsize,
int row)
359 for (
int j=0;
j<vertexsize;++
j) {
360 int indexj = set.subIndex(*it,
j,dim);
361 for (
int l=0;l<dim;++l) {
362 MPC* mpc = getMPC(indexj,l);
368 adjacencyPattern[row].insert(idx);
370 }
else if (meqn[indexj*dim+l] != -1)
371 adjacencyPattern[row].insert(meqn[indexj*dim+l]);
376 template<
class Gr
idType>
379 adjacencyPattern.resize(maxeqn);
380 std::cout <<
"\tsetting up sparsity pattern..." << std::endl;
388 for (
LeafIterator it = gv.leafGridView().template begin<0>(); it !=
itend; ++it, ++cell) {
389 Dune::GeometryType gt = it->type();
391 const Dune::template ReferenceElement<double,dim>& ref =
392 Dune::ReferenceElements<double,dim>::general(gt);
394 int vertexsize = ref.size(dim);
395 for (
int i=0; i < vertexsize; i++) {
396 int indexi = set.subIndex(*it,i,dim);
397 for (
int k=0;k<dim;++k) {
398 MPC* mpc = getMPC(indexi,k);
401 nodeAdjacency(it,vertexsize,
406 nodeAdjacency(it,vertexsize,meqn[indexi*dim+k]);
409 if (cell % 10000 == 0)
410 help.log(cell,
"\t\t... still processing ... cell ");
void removeMaster(size_t pos)
Removes the pos'th master DOF from the constraint equation.
Definition: mpc.hh:126
void printLoadVector() const
Print the current load vector.
Definition: asmhandler_impl.hpp:243
Definition: applier.hpp:18
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. Assemble entries for a single DOF.
Definition: asmhandler_impl.hpp:52
const DOF & getSlave() const
Returns a reference to the slave DOF.
Definition: mpc.hh:139
int node
Node number identifying this DOF.
Definition: mpc.hh:83
void preprocess()
Internal function. Generate meqn for registered MPC/fixed nodes.
Definition: asmhandler_impl.hpp:317
GridType::LeafGridView::template Codim< 0 >::Iterator LeafIterator
An iterator over grid cells.
Definition: asmhandler.hpp:38
void extractValues(Dune::FieldVector< double, comp > &v, const Vector &u, const LeafIterator &it)
Extract values corresponding to cell.
Definition: asmhandler_impl.hpp:114
A class for representing a general multi-point constraint equation.
Definition: mpc.hh:58
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
const DOF & getMaster(size_t i) const
Returns a reference to the i'th master DOF.
Definition: mpc.hh:141
void expandSolution(Vector &result, const Vector &u)
Expand a system vector to a solution vector.
Definition: asmhandler_impl.hpp:147
void resolveMPCChain(MPC *mpc)
Internal function. Handles a single MPC.
Definition: asmhandler_impl.hpp:252
size_t getNoMaster() const
Returns the number of master DOFs.
Definition: mpc.hh:144
static void fromAdjacency(Matrix &A, const AdjacencyPattern &adj, int rows, int cols)
Create a sparse matrix from a given adjacency pattern.
void updateFixedNode(int node, const std::pair< Direction, NodeValue > &entry)
Update/add a fixed node.
Definition: asmhandler_impl.hpp:216
void determineAdjacencyPattern()
Internal function. Calculate adjacency pattern.
Definition: asmhandler_impl.hpp:377
void printOperator() const
Print the current operator.
Definition: asmhandler_impl.hpp:237
void nodeAdjacency(const LeafIterator &it, int vertexsize, int row)
Internal function. Generate adjacency pattern for a given node.
Definition: asmhandler_impl.hpp:353
static void print(const Matrix &A)
Print a matrix to stdout.
void updateMaster(size_t pos, double c)
Updates the coefficient of the pos'th master DOF.
Definition: mpc.hh:119
int dof
Local DOF number within node.
Definition: mpc.hh:84
void addOffset(double offset)
Increments the c0 coefficient by a given offset.
Definition: mpc.hh:133
void initForAssembly()
This function needs to be called before starting the element assembly.
Definition: asmhandler_impl.hpp:22
Direction
An enum for specification of global coordinate directions.
Definition: mpc.hh:24
void addMPC(MPC *mpc)
Add a MPC.
Definition: asmhandler_impl.hpp:189
A struct for representing one term (DOF number and associated coefficient) in a MPC equation...
Definition: mpc.hh:66
double coeff
The constrained value, or master DOF scaling coefficient.
Definition: mpc.hh:85
MPC * getMPC(int node, int dof)
Look for and return a MPC for a specified node+dof.
Definition: asmhandler_impl.hpp:207
Dune::BlockVector< Dune::FieldVector< double, 1 > > Vector
A vector holding our RHS.
Definition: matrixops.hpp:29
const LeafVertexIterator itend
Definition: elasticity_upscale_impl.hpp:145
GridType::LeafGridView::IndexSet LeafIndexSet
A set of indices.
Definition: asmhandler.hpp:35
Helper class for progress logging during time consuming processes.
Definition: logutils.hpp:21
int j
Definition: elasticity_upscale_impl.hpp:658
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
fixMap::iterator fixIt
Iterator over a fixmap.
Definition: asmhandler.hpp:205