12#ifndef ASMHANDLER_IMPL_HPP_
13#define ASMHANDLER_IMPL_HPP_
15#include <dune/common/version.hh>
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 (*bptr)[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>
121 auto ref = Dune::ReferenceElements<double,dim>::cube();
122 int vertexsize = ref.size(dim);
124 for (
int i=0;i<vertexsize;++i) {
125 int indexi = set.subIndex(*it,i,dim);
126 fixIt it2 = fixedNodes.find(indexi);
127 for (
int n=0;n<dim;++n) {
128 MPC* mpc = getMPC(indexi,n);
129 if (it2 != fixedNodes.end() && it2->second.first & (1 << n))
130 v[l++] = it2->second.second[n];
139 v[l++] = u[meqn[indexi*dim+n]];
144 template<
class Gr
idType>
147 int nodes = gv.size(dim);
148 result.resize(nodes*dim);
151 for (
int i=0;i<nodes;++i) {
152 fixIt it = fixedNodes.find(i);
154 if (it == fixedNodes.end())
157 dir = it->second.first;
160 for (
int j=0;
j<dim;++
j) {
162 result[l] = it->second.second[
j];
163 else if (meqn[l] != -1)
164 result[l] = u[meqn[l]];
171 for (
int i=0;i<nodes;++i) {
172 for (
int j=0;
j<dim;++
j) {
173 MPC* mpc = getMPC(i,
j);
186 template<
class Gr
idType>
195 if (it == fixedNodes.end() ||
196 !(it->second.first & flag)) {
197 mpcs.insert(std::make_pair(slaveNode,mpc));
204 template<
class Gr
idType>
207 if (mpcs.find(node*dim+dof) != mpcs.end())
208 return mpcs[node*dim+dof];
213 template<
class Gr
idType>
215 const std::pair<Direction,NodeValue>& entry)
217 fixIt it = fixedNodes.find(node);
219 if (it == fixedNodes.end() || it->second.first == entry.first) {
220 fixedNodes[node] = entry;
223 int temp = it->second.first;
227 for (
int i=0;i<dim;++i) {
228 if (entry.first & flag)
229 it->second.second[i] = entry.second[i];
234 template<
class Gr
idType>
240 template<
class Gr
idType>
243 for (
int i=0;i<b.size();++i) {
245 std::cout << (fabs(val)>1.e-12?val:0.f) << std::endl;
249 template<
class Gr
idType>
253 if (nMaster == 0)
return;
255 for (
size_t i = 0; i < nMaster; i++)
260 Dune::FieldVector<double,dim> coeff;
265 for (
int d = 1; d <= dim; d++)
266 if (fabs(coeff[d-1]) > 1.0e-8)
274 resolveMPCChain(mpc2);
277 if (!removeOld) removeOld = 1;
314 template<
class Gr
idType>
317 int nodes = gv.size(dim);
318 meqn.resize(nodes*dim);
321 for (
int indexi=0;indexi<nodes;++indexi) {
322 fixIt it2 = fixedNodes.find(indexi);
323 if (it2 == fixedNodes.end()) {
324 for (
int i=0;i<dim;++i) {
325 MPC* mpc = getMPC(indexi,i);
327 meqn[indexi*dim+i] = maxeqn++;
329 meqn[indexi*dim+i] = -1;
333 for (
int i=0;i<dim;++i) {
334 if (it2->second.first & flag)
335 meqn[indexi*dim+i] = -1;
337 MPC* mpc = getMPC(indexi,i);
339 meqn[indexi*dim+i] = maxeqn++;
341 meqn[indexi*dim+i] = -1;
347 std::cout <<
"\tnumber of equations: " << maxeqn << std::endl;
350 template<
class Gr
idType>
352 int vertexsize,
int row)
357 for (
int j=0;
j<vertexsize;++
j) {
358 int indexj = set.subIndex(*it,
j,dim);
359 for (
int l=0;l<dim;++l) {
360 MPC* mpc = getMPC(indexj,l);
366 adjacencyPattern[row].insert(idx);
368 }
else if (meqn[indexj*dim+l] != -1)
369 adjacencyPattern[row].insert(meqn[indexj*dim+l]);
374 template<
class Gr
idType>
377 adjacencyPattern.resize(maxeqn);
378 std::cout <<
"\tsetting up sparsity pattern..." << std::endl;
386 for (
LeafIterator it = gv.leafGridView().template begin<0>(); it !=
itend; ++it, ++cell) {
387 auto ref = Dune::ReferenceElements<double,dim>::cube();
389 int vertexsize = ref.size(dim);
390 for (
int i=0; i < vertexsize; i++) {
391 int indexi = set.subIndex(*it,i,dim);
392 for (
int k=0;k<dim;++k) {
393 MPC* mpc = getMPC(indexi,k);
396 nodeAdjacency(it,vertexsize,
401 nodeAdjacency(it,vertexsize,meqn[indexi*dim+k]);
404 if (cell % 10000 == 0)
405 help.
log(cell,
"\t\t... still processing ... cell ");
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
void initForAssembly()
This function needs to be called before starting the element assembly.
Definition: asmhandler_impl.hpp:22
void resolveMPCChain(MPC *mpc)
Internal function. Handles a single MPC.
Definition: asmhandler_impl.hpp:250
fixMap::iterator fixIt
Iterator over a fixmap.
Definition: asmhandler.hpp:211
GridType::LeafGridView::IndexSet LeafIndexSet
A set of indices.
Definition: asmhandler.hpp:41
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
void updateFixedNode(int node, const std::pair< Direction, NodeValue > &entry)
Update/add a fixed node.
Definition: asmhandler_impl.hpp:214
void extractValues(Dune::FieldVector< double, comp > &v, const Vector &u, const LeafIterator &it)
Extract values corresponding to cell.
Definition: asmhandler_impl.hpp:114
void expandSolution(Vector &result, const Vector &u)
Expand a system vector to a solution vector.
Definition: asmhandler_impl.hpp:145
void addMPC(MPC *mpc)
Add a MPC.
Definition: asmhandler_impl.hpp:187
void preprocess()
Internal function. Generate meqn for registered MPC/fixed nodes.
Definition: asmhandler_impl.hpp:315
void nodeAdjacency(const LeafIterator &it, int vertexsize, int row)
Internal function. Generate adjacency pattern for a given node.
Definition: asmhandler_impl.hpp:351
GridType::LeafGridView::template Codim< 0 >::Iterator LeafIterator
An iterator over grid cells.
Definition: asmhandler.hpp:44
void printLoadVector() const
Print the current load vector.
Definition: asmhandler_impl.hpp:241
MPC * getMPC(int node, int dof)
Look for and return a MPC for a specified node+dof.
Definition: asmhandler_impl.hpp:205
void determineAdjacencyPattern()
Internal function. Calculate adjacency pattern.
Definition: asmhandler_impl.hpp:375
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
void printOperator() const
Print the current operator.
Definition: asmhandler_impl.hpp:235
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 & getSlave() const
Returns a reference to the slave DOF.
Definition: mpc.hh:136
void addOffset(double offset)
Increments the c0 coefficient by a given offset.
Definition: mpc.hh:130
void removeMaster(size_t pos)
Removes the pos'th master DOF from the constraint equation.
Definition: mpc.hh:123
void updateMaster(size_t pos, double c)
Updates the coefficient of the pos'th master DOF.
Definition: mpc.hh:116
const DOF & getMaster(size_t i) const
Returns a reference to the i'th master DOF.
Definition: mpc.hh:138
static void fromAdjacency(Matrix &A, const AdjacencyPattern &adj, int rows, int cols)
Create a sparse matrix from a given adjacency pattern.
static void print(const Matrix &A)
Print a matrix to stdout.
Direction
An enum for specification of global coordinate directions.
Definition: mpc.hh:24
@ XYZ
Definition: mpc.hh:26
@ NONE
Definition: mpc.hh:24
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
Definition: ImplicitAssembly.hpp:43
A struct for representing one term (DOF number and associated coefficient) in a MPC equation.
Definition: mpc.hh:67
int dof
Local DOF number within node.
Definition: mpc.hh:81
int node
Node number identifying this DOF.
Definition: mpc.hh:80
double coeff
The constrained value, or master DOF scaling coefficient.
Definition: mpc.hh:82