My Project
asmhandler_impl.hpp
Go to the documentation of this file.
1 //==============================================================================
11 //==============================================================================
12 #ifndef ASMHANDLER_IMPL_HPP_
13 #define ASMHANDLER_IMPL_HPP_
14 
15 #include <dune/common/version.hh>
16 #include <iostream>
17 
18 namespace Opm {
19 namespace Elasticity {
20 
21  template<class GridType>
23 {
24  resolveMPCChains();
25  preprocess();
26  determineAdjacencyPattern();
27 
28  MatrixOps::fromAdjacency(A,adjacencyPattern,
29  adjacencyPattern.size(),adjacencyPattern.size());
30  b.resize(adjacencyPattern.size());
31  b = 0;
32  adjacencyPattern.clear();
33 
34  // print some information
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;
38  int fixedDofs=0;
39  for (fixIt it = fixedNodes.begin(); it != fixedNodes.end(); ++it) {
40  if (it->second.first & X)
41  fixedDofs++;
42  if (it->second.first & Y)
43  fixedDofs++;
44  if (it->second.first & Z)
45  fixedDofs++;
46  }
47  std::cout << "\tNumber of fixed dofs: " << fixedDofs << std::endl;
48 }
49 
50  template<class GridType>
51  template<int esize>
52 void ASMHandler<GridType>::addDOF(int row, int erow,
53  const Dune::FieldMatrix<double,esize,esize>* K,
54  const Dune::FieldVector<double,esize>* S,
55  const LeafIndexSet& set,
56  const LeafIterator& cell,
57  Vector* bptr,
58  double scale)
59 {
60  if (row == -1)
61  return;
62  if (K) {
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);
67  if (mpc) {
68  for (size_t n=0;n<mpc->getNoMaster();++n) {
69  int idx = meqn[mpc->getMaster(n).node*dim+mpc->getMaster(n).dof-1];
70  if (idx != -1)
71  A[row][idx] += scale*mpc->getMaster(n).coeff*(*K)[erow][j*dim+l];
72  }
73  } else if (meqn[index2*dim+l] != -1) {
74  A[row][meqn[index2*dim+l]] += scale*(*K)[erow][j*dim+l];
75  }
76  }
77  }
78  }
79  if (S && bptr)
80  (*bptr)[row] += scale*(*S)[erow];
81 }
82 
83  template<class GridType>
84  template<int esize>
86  const Dune::FieldMatrix<double,esize,esize>* K,
87  const Dune::FieldVector<double,esize>* S,
88  const LeafIterator& cell,
89  Vector* b2)
90 {
91  if (!b2)
92  b2 = &b;
93  const LeafIndexSet& set = gv.leafGridView().indexSet();
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)
98  continue;
99  for (int k=0;k<dim;++k) {
100  MPC* mpc = getMPC(index1,k);
101  if (mpc) {
102  for (size_t n=0;n<mpc->getNoMaster();++n) {
103  int idx = meqn[mpc->getMaster(n).node*dim+mpc->getMaster(n).dof-1];
104  addDOF(idx,i*dim+k,K,S,set,cell,b2,mpc->getMaster(n).coeff);
105  }
106  } else
107  addDOF(meqn[index1*dim+k],i*dim+k,K,S,set,cell,b2);
108  }
109  }
110 }
111 
112  template<class GridType>
113  template<int comp>
114 void ASMHandler<GridType>::extractValues(Dune::FieldVector<double,comp>& v,
115  const Vector& u,
116  const LeafIterator& it)
117 {
118  v = 0;
119  const LeafIndexSet& set = gv.leafGridView().indexSet();
120 
121  auto ref = Dune::ReferenceElements<double,dim>::cube();
122  int vertexsize = ref.size(dim);
123  int l=0;
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];
131  else if (mpc) {
132  for (size_t m=0;m<mpc->getNoMaster();++m) {
133  int idx = meqn[mpc->getMaster(m).node*dim+mpc->getMaster(m).dof-1];
134  if (idx != -1)
135  v[l] += u[idx]*mpc->getMaster(m).coeff;
136  }
137  l++;
138  } else
139  v[l++] = u[meqn[indexi*dim+n]];
140  }
141  }
142 }
143 
144  template<class GridType>
146 {
147  int nodes = gv.size(dim);
148  result.resize(nodes*dim);
149  result = 0;
150  int l=0;
151  for (int i=0;i<nodes;++i) {
152  fixIt it = fixedNodes.find(i);
153  Direction dir;
154  if (it == fixedNodes.end())
155  dir = NONE;
156  else
157  dir = it->second.first;
158 
159  int flag=1;
160  for (int j=0;j<dim;++j) {
161  if (dir & flag)
162  result[l] = it->second.second[j];
163  else if (meqn[l] != -1)
164  result[l] = u[meqn[l]];
165  l++;
166  flag *= 2;
167  }
168  }
169  // second loop - handle MPC couplings
170  l = 0;
171  for (int i=0;i<nodes;++i) {
172  for (int j=0;j<dim;++j) {
173  MPC* mpc = getMPC(i,j);
174  if (mpc) {
175  for (size_t n=0;n<mpc->getNoMaster();++n) {
176  int idx = mpc->getMaster(n).node*dim+mpc->getMaster(n).dof-1;
177  if (meqn[idx] != -1)
178  result[l] += u[meqn[idx]]*mpc->getMaster(n).coeff;
179  }
180  }
181  ++l;
182  }
183  }
184 }
185 
186  template<class GridType>
188 {
189  if (!mpc)
190  return;
191 
192  int slaveNode = mpc->getSlave().node*dim+mpc->getSlave().dof-1;
193  fixIt it = fixedNodes.find(mpc->getSlave().node);
194  int flag = 1 << (mpc->getSlave().dof-1);
195  if (it == fixedNodes.end() ||
196  !(it->second.first & flag)) {
197  mpcs.insert(std::make_pair(slaveNode,mpc));
198  return;
199  }
200 
201  delete mpc;
202 }
203 
204  template<class GridType>
206 {
207  if (mpcs.find(node*dim+dof) != mpcs.end())
208  return mpcs[node*dim+dof];
209 
210  return NULL;
211 }
212 
213  template<class GridType>
215  const std::pair<Direction,NodeValue>& entry)
216 {
217  fixIt it = fixedNodes.find(node);
218  // same type or new - update values/add and return
219  if (it == fixedNodes.end() || it->second.first == entry.first) {
220  fixedNodes[node] = entry;
221  return;
222  }
223  int temp = it->second.first;
224  temp |= entry.first;
225  it->second.first = (Direction)temp;
226  int flag = 1;
227  for (int i=0;i<dim;++i) {
228  if (entry.first & flag)
229  it->second.second[i] = entry.second[i];
230  flag *= 2;
231  }
232 }
233 
234  template<class GridType>
236 {
237  MatrixOps::print(A);
238 }
239 
240  template<class GridType>
242 {
243  for (int i=0;i<b.size();++i) {
244  double val = b[i];
245  std::cout << (fabs(val)>1.e-12?val:0.f) << std::endl;
246  }
247 }
248 
249  template<class GridType>
251 {
252  size_t nMaster = mpc->getNoMaster();
253  if (nMaster == 0) return; // no masters, prescribed displacement only
254 
255  for (size_t i = 0; i < nMaster; i++)
256  {
257  // Check if the master node has a local coordinate system attached. If yes,
258  // the slave DOF might be coupled to all (local) DOFs of the master node.
259  const MPC::DOF& master = mpc->getMaster(i);
260  Dune::FieldVector<double,dim> coeff;
261  coeff = 0;
262  coeff[master.dof-1] = master.coeff;
263 
264  int removeOld = 0;
265  for (int d = 1; d <= dim; d++)
266  if (fabs(coeff[d-1]) > 1.0e-8)
267  {
268  MPC* mpc2 = getMPC(mpc->getMaster(i).node,d-1);
269  if (mpc2)
270  {
271  // We have a master DOF which is a slave in another MPC.
272  // Invoke resolveMPCchain recursively to ensure that all master DOFs
273  // of that equation are not slaves themselves.
274  resolveMPCChain(mpc2);
275 
276  // Remove current master specification, unless it has been updated
277  if (!removeOld) removeOld = 1;
278 
279  // Add constant offset from the other equation
280  mpc->addOffset(mpc2->getSlave().coeff*coeff[d-1]);
281 
282  // Add masters from the other equations
283  for (size_t j = 0; j < mpc2->getNoMaster(); j++)
284  mpc->addMaster(mpc2->getMaster(j).node,
285  mpc2->getMaster(j).dof,
286  mpc2->getMaster(j).coeff*coeff[d-1]);
287  }
288  else
289  // The master node is free, but has a local coordinate system
290  if (d != mpc->getMaster(i).dof)
291  // Add coupling to the other local DOFs of this master node.
292  mpc->addMaster(mpc->getMaster(i).node,d,coeff[d-1]);
293  else if (coeff[d-1] != mpc->getMaster(i).coeff)
294  {
295  // Update the coupling coefficient of this master DOF
296  // due to the local-to-global transformation
297  mpc->updateMaster(i,coeff[d-1]);
298  removeOld = -1;
299  }
300  }
301  else if (d == mpc->getMaster(i).dof && !removeOld)
302  // The coefficient of the current master DOF is zero,
303  removeOld = 1; // so remove it from the contraint equation
304 
305  if (removeOld == 1)
306  {
307  // Remove the old master DOF specification when it has been replaced
308  mpc->removeMaster(i--);
309  nMaster--; // we don't need to check the added masters
310  }
311  }
312 }
313 
314  template<class GridType>
316 {
317  int nodes = gv.size(dim);
318  meqn.resize(nodes*dim);
319 
320  // iterate over nodes
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);
326  if (!mpc)
327  meqn[indexi*dim+i] = maxeqn++;
328  else
329  meqn[indexi*dim+i] = -1;
330  }
331  } else {
332  int flag=1;
333  for (int i=0;i<dim;++i) {
334  if (it2->second.first & flag)
335  meqn[indexi*dim+i] = -1;
336  else {
337  MPC* mpc = getMPC(indexi,i);
338  if (!mpc)
339  meqn[indexi*dim+i] = maxeqn++;
340  else
341  meqn[indexi*dim+i] = -1;
342  }
343  flag *= 2;
344  }
345  }
346  }
347  std::cout << "\tnumber of equations: " << maxeqn << std::endl;
348 }
349 
350  template<class GridType>
352  int vertexsize, int row)
353 {
354  if (row == -1)
355  return;
356  const LeafIndexSet& set = gv.leafGridView().indexSet();
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);
361  if (mpc) {
362  for (size_t i=0;i<mpc->getNoMaster();++i) {
363  int idx = meqn[mpc->getMaster(i).node*dim+
364  mpc->getMaster(i).dof-1];
365  if (idx != -1)
366  adjacencyPattern[row].insert(idx);
367  }
368  } else if (meqn[indexj*dim+l] != -1)
369  adjacencyPattern[row].insert(meqn[indexj*dim+l]);
370  }
371  }
372 }
373 
374  template<class GridType>
376 {
377  adjacencyPattern.resize(maxeqn);
378  std::cout << "\tsetting up sparsity pattern..." << std::endl;
379  LoggerHelper help(gv.size(0), 5, 50000);
380 
381  const LeafIndexSet& set = gv.leafGridView().indexSet();
382  LeafIterator itend = gv.leafGridView().template end<0>();
383 
384  // iterate over cells
385  int cell=0;
386  for (LeafIterator it = gv.leafGridView().template begin<0>(); it != itend; ++it, ++cell) {
387  auto ref = Dune::ReferenceElements<double,dim>::cube();
388 
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);
394  if (mpc) {
395  for (size_t l=0;l<mpc->getNoMaster();++l) {
396  nodeAdjacency(it,vertexsize,
397  meqn[mpc->getMaster(l).node*dim+
398  mpc->getMaster(l).dof-1]);
399  }
400  } else
401  nodeAdjacency(it,vertexsize,meqn[indexi*dim+k]);
402  }
403  }
404  if (cell % 10000 == 0)
405  help.log(cell, "\t\t... still processing ... cell ");
406  }
407 }
408 
409 }} // namespace Opm, Elasticity
410 
411 #endif
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.
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.
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.
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:144
const DOF & getMaster(size_t i) const
Returns a reference to the i'th master DOF.
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:113
const DOF & getSlave() const
Returns a reference to the slave DOF.
Definition: mpc.hh:139
void addOffset(double offset)
Increments the c0 coefficient by a given offset.
Definition: mpc.hh:133
void removeMaster(size_t pos)
Removes the pos'th master DOF from the constraint equation.
Definition: mpc.hh:126
void updateMaster(size_t pos, double c)
Updates the coefficient of the pos'th master DOF.
Definition: mpc.hh:119
static void print(const Matrix &A)
Print a matrix to stdout.
Definition: matrixops.cpp:73
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
Dune::BlockVector< Dune::FieldVector< double, 1 > > Vector
A vector holding our RHS.
Definition: matrixops.hpp:33
Direction
An enum for specification of global coordinate directions.
Definition: mpc.hh:24
Inverting small matrices.
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:84
int node
Node number identifying this DOF.
Definition: mpc.hh:83
double coeff
The constrained value, or master DOF scaling coefficient.
Definition: mpc.hh:85