1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 
6 #ifndef DOF_MANAGER_H
7 #define DOF_MANAGER_H
8 
9 #include <vector>
10 #include <string>
11 #include <complex>
12 #include <map>
13 #include <list>
14 #include <iostream>
15 #include "MVertex.h"
16 #include "linearSystem.h"
17 #include "fullMatrix.h"
18 
19 class Dof {
20 protected:
21   // v(x) = \sum_f \sum_i v_{fi} s^f_i(x)
22   long int _entity; // "i": node, edge, group, etc.
23   int _type; // "f": basis function type index, etc.
24 public:
Dof(long int entity,int type)25   Dof(long int entity, int type) : _entity(entity), _type(type) {}
getEntity()26   inline long int getEntity() const { return _entity; }
getType()27   inline int getType() const { return _type; }
createTypeWithTwoInts(int i1,int i2)28   inline static int createTypeWithTwoInts(int i1, int i2)
29   {
30     return i1 + 10000 * i2;
31   }
getTwoIntsFromType(int t,int & i1,int & i2)32   inline static void getTwoIntsFromType(int t, int &i1, int &i2)
33   {
34     i1 = t % 10000;
35     i2 = t / 10000;
36   }
37   bool operator<(const Dof &other) const
38   {
39     if(_entity < other._entity) return true;
40     if(_entity > other._entity) return false;
41     if(_type < other._type) return true;
42     return false;
43   }
44   bool operator==(const Dof &other) const
45   {
46     return (_entity == other._entity && _type == other._type);
47   }
48 };
49 
50 template <class T> struct dofTraits {
51   typedef T VecType;
52   typedef T MatType;
gemmdofTraits53   inline static void gemm(VecType &r, const MatType &m, const VecType &v,
54                           double alpha, double beta)
55   {
56     r = beta * r + alpha * m * v;
57   }
58 };
59 
60 template <class T> struct dofTraits<fullMatrix<T> > {
61   typedef fullMatrix<T> VecType;
62   typedef fullMatrix<T> MatType;
63   inline static void gemm(VecType &r, const MatType &m, const VecType &v,
64                           double alpha, double beta)
65   {
66     r.gemm(m, v, alpha, beta);
67   }
68 };
69 
70 /*
71 template<> struct dofTraits<fullVector<std::complex<double> > >
72 {
73   typedef fullVector<std::complex<double> > VecType;
74   typedef fullMatrix<std::complex<double> > MatType;
75 };
76 */
77 
78 template <class T> class DofAffineConstraint {
79 public:
80   std::vector<std::pair<Dof, typename dofTraits<T>::MatType> > linear;
81   typename dofTraits<T>::VecType shift;
82 };
83 
84 // non template part that can be implemented in the cxx file (and so avoid to
85 // include mpi.h in the .h file)
86 class dofManagerBase {
87 protected:
88   // numbering of unknown dof blocks
89   std::map<Dof, int> unknown;
90 
91   // associatations (not used ?)
92   std::map<Dof, Dof> associatedWith;
93 
94   // parallel section
95   // those dof are images of ghost located on another proc (id givent by the
96   // map). this is a first try, maybe not the final implementation
97   std::map<Dof, std::pair<int, int> > ghostByDof; // dof => procId, globalId
98   std::vector<std::vector<Dof> > ghostByProc, parentByProc;
99   int _localSize;
100   bool _parallelFinalized;
101   bool _isParallel;
102   void _parallelFinalize();
103   dofManagerBase(bool isParallel)
104   {
105     _isParallel = isParallel;
106     _parallelFinalized = false;
107   }
108 };
109 
110 // A manager for degrees of freedoms, templated on the value of a dof
111 // (what the functional returns): float, double, complex<double>,
112 // fullVecor<double>, ...
113 template <class T> class dofManager : public dofManagerBase {
114 public:
115   typedef typename dofTraits<T>::VecType dataVec;
116   typedef typename dofTraits<T>::MatType dataMat;
117 
118 protected:
119   // general affine constraint on sub-blocks, treated by adding
120   // equations:
121   //   Dof = \sum_i dataMat_i x Dof_i + dataVec
122   std::map<Dof, DofAffineConstraint<dataVec> > constraints;
123 
124   // fixations on full blocks, treated by eliminating equations:
125   //   DofVec = dataVec
126   std::map<Dof, dataVec> fixed;
127 
128   // initial conditions (not used ?)
129   std::map<Dof, std::vector<dataVec> > initial;
130 
131   // linearSystems
132   linearSystem<dataMat> *_current;
133   std::map<const std::string, linearSystem<dataMat> *> _linearSystems;
134 
135   std::map<Dof, T> ghostValue;
136 
137 public:
138   void scatterSolution();
139 
140 public:
141   dofManager(linearSystem<dataMat> *l, bool isParallel = false)
142     : dofManagerBase(isParallel), _current(l)
143   {
144     _linearSystems["A"] = l;
145   }
146   dofManager(linearSystem<dataMat> *l1, linearSystem<dataMat> *l2)
147     : dofManagerBase(false), _current(l1)
148   {
149     _linearSystems.insert(std::make_pair("A", l1));
150     _linearSystems.insert(std::make_pair("B", l2));
151   }
152   virtual ~dofManager() {}
153   virtual inline void fixDof(Dof key, const dataVec &value)
154   {
155     if(unknown.find(key) != unknown.end()) return;
156     fixed[key] = value;
157   }
158   inline void fixDof(long int ent, int type, const dataVec &value)
159   {
160     fixDof(Dof(ent, type), value);
161   }
162   inline void associateDof(long int ent_from, int type_from,
163 			   long int ent_to, int type_to)
164   {
165     Dof from (ent_from, type_from);
166     Dof to   (ent_to, type_to);
167     std::pair<Dof,Dof> p = std::make_pair(from,to);
168     associatedWith.insert(p);
169   }
170   void fixVertex(MVertex const *v, int iComp, int iField, const dataVec &value)
171   {
172     fixDof(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField), value);
173   }
174   virtual inline bool isFixed(Dof key) const
175   {
176     if(fixed.find(key) != fixed.end()) {
177       return true;
178     }
179     return false;
180   }
181 
182   virtual inline bool isAnUnknown(Dof key) const
183   {
184     if(ghostValue.find(key) == ghostValue.end()) {
185       if(unknown.find(key) != unknown.end()) return true;
186     }
187     return false;
188   }
189 
190   virtual inline bool isConstrained(Dof key) const
191   {
192     if(constraints.find(key) != constraints.end()) {
193       return true;
194     }
195     return false;
196   }
197 
198   inline bool isFixed(long int ent, int type) const
199   {
200     return isFixed(Dof(ent, type));
201   }
202   inline bool isFixed(MVertex *v, int iComp, int iField) const
203   {
204     return isFixed(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField));
205   }
206   virtual inline void numberGhostDof(Dof key, int procId)
207   {
208     if(fixed.find(key) != fixed.end()) return;
209     if(constraints.find(key) != constraints.end()) return;
210     if(ghostByDof.find(key) != ghostByDof.end()) return;
211     ghostByDof[key] = std::make_pair(procId, 0);
212   }
213   virtual inline void numberDof(Dof key)
214   {
215     if(associatedWith.find(key) != associatedWith.end()) return;
216     if(fixed.find(key) != fixed.end()) return;
217     if(constraints.find(key) != constraints.end()) return;
218     if(ghostByDof.find(key) != ghostByDof.end()) return;
219 
220     auto it = unknown.find(key);
221     if(it == unknown.end()) {
222       std::size_t size = unknown.size();
223       unknown[key] = size;
224     }
225   }
226   virtual inline void numberDof(const std::vector<Dof> &R)
227   {
228     for(std::size_t i = 0; i < R.size(); i++) this->numberDof(R[i]);
229   }
230   inline void numberDof(long int ent, int type) { numberDof(Dof(ent, type)); }
231   inline void numberVertex(MVertex *v, int iComp, int iField)
232   {
233     numberDof(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField));
234   }
235   virtual inline void getDofValue(std::vector<Dof> &keys,
236                                   std::vector<dataVec> &Vals)
237   {
238     for(std::size_t i = 0; i < keys.size(); i++) {
239       auto it = associatedWith.find(keys[i]);
240       if (it != associatedWith.end())keys[i] = it->second;
241     }
242 
243     int ndofs = keys.size();
244     size_t originalSize = Vals.size();
245     Vals.resize(originalSize + ndofs);
246     for(int i = 0; i < ndofs; ++i) getDofValue(keys[i], Vals[originalSize + i]);
247   }
248 
249   virtual inline bool getAnUnknown(Dof key, dataVec &val) const
250   {
251     if(ghostValue.find(key) == ghostValue.end()) {
252       auto it = unknown.find(key);
253       if(it != unknown.end()) {
254         _current->getFromSolution(it->second, val);
255         return true;
256       }
257     }
258     return false;
259   }
260 
261   virtual inline void getFixedDofValue(Dof key, dataVec &val) const
262   {
263     typename std::map<Dof, dataVec>::const_iterator it = fixed.find(key);
264     if(it != fixed.end()) {
265       val = it->second;
266     }
267     else {
268       Msg::Error("getFixedDof: Dof is not fixed");
269       return;
270     }
271   };
272 
273   virtual inline void getDofValue(Dof key, dataVec &val) const
274   {
275     {
276       auto it = associatedWith.find(key);
277       if (it != associatedWith.end()){
278 	//	  printf("ass to %d\n",it->second.getEntity());
279 	auto itx = unknown.find(it->second);
280 	if(itx != unknown.end()) {
281 	  _current->getFromSolution(itx->second, val);
282 	  return;
283 	}
284 	key = it->second;
285       }
286     }
287     {
288       typename std::map<Dof, dataVec>::const_iterator it = ghostValue.find(key);
289       if(it != ghostValue.end()) {
290         val = it->second;
291         return;
292       }
293     }
294     {
295       auto it = unknown.find(key);
296       if(it != unknown.end()) {
297         _current->getFromSolution(it->second, val);
298         return;
299       }
300     }
301     {
302       typename std::map<Dof, dataVec>::const_iterator it = fixed.find(key);
303       if(it != fixed.end()) {
304         val = it->second;
305         return;
306       }
307     }
308     {
309       typename std::map<Dof, DofAffineConstraint<dataVec> >::const_iterator it =
310         constraints.find(key);
311       if(it != constraints.end()) {
312         dataVec tmp(val);
313         val = it->second.shift;
314         for(unsigned i = 0; i < (it->second).linear.size(); i++) {
315           /* gcc: warning: variable ‘itu’ set but not used
316           std::map<Dof, int>::const_iterator itu = unknown.find
317             (((it->second).linear[i]).first);*/
318           getDofValue(((it->second).linear[i]).first, tmp);
319           dofTraits<T>::gemm(val, ((it->second).linear[i]).second, tmp, 1, 1);
320         }
321         return;
322       }
323     }
324   }
325   inline void getDofValue(int ent, int type, dataVec &v) const
326   {
327     getDofValue(Dof(ent, type), v);
328   }
329   inline void getDofValue(MVertex *v, int iComp, int iField,
330                           dataVec &value) const
331   {
332     getDofValue(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField), value);
333   }
334 
335   virtual inline void insertInSparsityPatternLinConst(const Dof &R,
336                                                       const Dof &C)
337   {
338     auto itR = unknown.find(R);
339     if(itR != unknown.end()) {
340       typename std::map<Dof, DofAffineConstraint<dataVec> >::iterator
341         itConstraint;
342       itConstraint = constraints.find(C);
343       if(itConstraint != constraints.end()) {
344         for(unsigned i = 0; i < (itConstraint->second).linear.size(); i++) {
345           insertInSparsityPattern(R, (itConstraint->second).linear[i].first);
346         }
347       }
348     }
349     else { // test function ; (no shift ?)
350       typename std::map<Dof, DofAffineConstraint<dataVec> >::iterator
351         itConstraint;
352       itConstraint = constraints.find(R);
353       if(itConstraint != constraints.end()) {
354         for(unsigned i = 0; i < (itConstraint->second).linear.size(); i++) {
355           insertInSparsityPattern((itConstraint->second).linear[i].first, C);
356         }
357       }
358     }
359   }
360 
361   virtual inline void insertInSparsityPattern(const Dof &R, const Dof &C)
362   {
363     if(_isParallel && !_parallelFinalized) _parallelFinalize();
364     if(!_current->isAllocated()) _current->allocate(sizeOfR());
365     auto itR = unknown.find(R);
366     if(itR != unknown.end()) {
367       auto itC = unknown.find(C);
368       if(itC != unknown.end()) {
369         _current->insertInSparsityPattern(itR->second, itC->second);
370       }
371       else {
372         typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(C);
373         if(itFixed != fixed.end()) {
374         }
375         else
376           insertInSparsityPatternLinConst(R, C);
377       }
378     }
379     if(itR == unknown.end()) {
380       insertInSparsityPatternLinConst(R, C);
381     }
382   }
383 
384   virtual inline void sparsityDof(const std::vector<Dof> &keys)
385   {
386     for(std::size_t itR = 0; itR < keys.size(); itR++) {
387       for(std::size_t itC = 0; itC < keys.size(); itC++) {
388         insertInSparsityPattern(keys[itR], keys[itC]);
389       }
390     }
391   }
392 
393   virtual inline void assemble(const Dof &R, const Dof &C, const dataMat &value)
394   {
395     if(_isParallel && !_parallelFinalized) _parallelFinalize();
396     if(!_current->isAllocated()) _current->allocate(sizeOfR());
397     auto itR = unknown.find(R);
398     if(itR != unknown.end()) {
399       auto itC = unknown.find(C);
400       if(itC != unknown.end()) {
401         _current->addToMatrix(itR->second, itC->second, value);
402       }
403       else {
404         typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(C);
405         if(itFixed != fixed.end()) {
406           // tmp = -value * itFixed->second
407           dataVec tmp(itFixed->second);
408           dofTraits<T>::gemm(tmp, value, itFixed->second, -1, 0);
409           _current->addToRightHandSide(itR->second, tmp);
410         }
411         else
412           assembleLinConst(R, C, value);
413       }
414     }
415     if(itR == unknown.end()) {
416       assembleLinConst(R, C, value);
417     }
418   }
419   virtual inline void assemble(std::vector<Dof> &R, std::vector<Dof> &C,
420                                const fullMatrix<dataMat> &m)
421   {
422     if(_isParallel && !_parallelFinalized) _parallelFinalize();
423     if(!_current->isAllocated()) _current->allocate(sizeOfR());
424     printf("coucou\n");
425 
426     for(std::size_t i = 0; i < R.size(); i++) {
427       auto it = associatedWith.find(R[i]);
428       if (it != associatedWith.end())R[i] = it->second;
429     }
430     for(std::size_t i = 0; i < C.size(); i++) {
431       auto it = associatedWith.find(C[i]);
432       if (it != associatedWith.end())C[i] = it->second;
433     }
434 
435     std::vector<int> NR(R.size()), NC(C.size());
436 
437     for(std::size_t i = 0; i < R.size(); i++) {
438       auto itR = unknown.find(R[i]);
439       if(itR != unknown.end())
440         NR[i] = itR->second;
441       else
442         NR[i] = -1;
443     }
444     for(std::size_t i = 0; i < C.size(); i++) {
445       auto itC = unknown.find(C[i]);
446       if(itC != unknown.end())
447         NC[i] = itC->second;
448       else
449         NC[i] = -1;
450     }
451     for(std::size_t i = 0; i < R.size(); i++) {
452       if(NR[i] != -1) {
453         for(std::size_t j = 0; j < C.size(); j++) {
454           if(NC[j] != -1) {
455             _current->addToMatrix(NR[i], NC[j], m(i, j));
456           }
457           else {
458             typename std::map<Dof, dataVec>::iterator itFixed =
459               fixed.find(C[j]);
460             if(itFixed != fixed.end()) {
461               // tmp = -m(i,j) * itFixed->second
462               dataVec tmp(itFixed->second);
463               dofTraits<T>::gemm(tmp, m(i, j), itFixed->second, -1, 0);
464               _current->addToRightHandSide(NR[i], tmp);
465             }
466             else
467               assembleLinConst(R[i], C[j], m(i, j));
468           }
469         }
470       }
471       else {
472         for(std::size_t j = 0; j < C.size(); j++) {
473           assembleLinConst(R[i], C[j], m(i, j));
474         }
475       }
476     }
477   }
478   // for linear forms
479   virtual inline void assemble(std::vector<Dof> &R,
480                                const fullVector<dataMat> &m)
481   {
482     if(_isParallel && !_parallelFinalized) _parallelFinalize();
483     if(!_current->isAllocated()) _current->allocate(sizeOfR());
484     printf("coucou RHS\n");
485 
486     for(std::size_t i = 0; i < R.size(); i++) {
487       auto it = associatedWith.find(R[i]);
488       if (it != associatedWith.end())R[i] = it->second;
489     }
490 
491 
492     std::vector<int> NR(R.size());
493     for(std::size_t i = 0; i < R.size(); i++) {
494       auto itR = unknown.find(R[i]);
495       if(itR != unknown.end())
496         NR[i] = itR->second;
497       else
498         NR[i] = -1;
499     }
500     for(std::size_t i = 0; i < R.size(); i++) {
501       if(NR[i] != -1) {
502         _current->addToRightHandSide(NR[i], m(i));
503       }
504       else {
505         typename std::map<Dof, DofAffineConstraint<dataVec> >::iterator
506           itConstraint;
507         itConstraint = constraints.find(R[i]);
508         if(itConstraint != constraints.end()) {
509           for(unsigned j = 0; j < (itConstraint->second).linear.size(); j++) {
510             dataMat tmp;
511             dofTraits<T>::gemm(tmp, (itConstraint->second).linear[j].second,
512                                m(i), 1, 0);
513             assemble((itConstraint->second).linear[j].first, tmp);
514           }
515         }
516       }
517     }
518   }
519   virtual inline void assemble(std::vector<Dof> &R,
520                                const fullMatrix<dataMat> &m)
521   {
522     if(_isParallel && !_parallelFinalized) _parallelFinalize();
523     if(!_current->isAllocated()) _current->allocate(sizeOfR());
524     for(std::size_t i = 0; i < R.size(); i++) {
525       auto it = associatedWith.find(R[i]);
526       if (it != associatedWith.end())R[i] = it->second;
527     }
528 
529     std::vector<int> NR(R.size());
530     for(std::size_t i = 0; i < R.size(); i++) {
531       auto itR = unknown.find(R[i]);
532       if(itR != unknown.end())
533         NR[i] = itR->second;
534       else
535         NR[i] = -1;
536     }
537     for(std::size_t i = 0; i < R.size(); i++) {
538       if(NR[i] != -1) {
539         for(std::size_t j = 0; j < R.size(); j++) {
540           if(NR[j] != -1) {
541             _current->addToMatrix(NR[i], NR[j], m(i, j));
542           }
543           else {
544             typename std::map<Dof, dataVec>::iterator itFixed =
545               fixed.find(R[j]);
546             if(itFixed != fixed.end()) {
547               // tmp = -m(i,j) * itFixed->second
548               dataVec tmp(itFixed->second);
549               dofTraits<T>::gemm(tmp, m(i, j), itFixed->second, -1, 0);
550               _current->addToRightHandSide(NR[i], tmp);
551             }
552             else
553               assembleLinConst(R[i], R[j], m(i, j));
554           }
555         }
556       }
557       else {
558         for(std::size_t j = 0; j < R.size(); j++) {
559           assembleLinConst(R[i], R[j], m(i, j));
560         }
561       }
562     }
563   }
564   inline void assemble(int entR, int typeR, int entC, int typeC,
565                        const dataMat &value)
566   {
567     assemble(Dof(entR, typeR), Dof(entC, typeC), value);
568   }
569   inline void assemble(MVertex *vR, int iCompR, int iFieldR, MVertex *vC,
570                        int iCompC, int iFieldC, const dataMat &value)
571   {
572     assemble(vR->getNum(), Dof::createTypeWithTwoInts(iCompR, iFieldR),
573              vC->getNum(), Dof::createTypeWithTwoInts(iCompC, iFieldC), value);
574   }
575   virtual inline void assemble(const Dof &R, const dataMat &value)
576   {
577     if(_isParallel && !_parallelFinalized) _parallelFinalize();
578     if(!_current->isAllocated()) _current->allocate(sizeOfR());
579     auto itR = unknown.find(R);
580     if(itR != unknown.end()) {
581       _current->addToRightHandSide(itR->second, value);
582     }
583     else {
584       typename std::map<Dof, DofAffineConstraint<dataVec> >::iterator
585         itConstraint;
586       itConstraint = constraints.find(R);
587       if(itConstraint != constraints.end()) {
588         for(unsigned j = 0; j < (itConstraint->second).linear.size(); j++) {
589           dataMat tmp;
590           dofTraits<T>::gemm(tmp, (itConstraint->second).linear[j].second,
591                              value, 1, 0);
592           assemble((itConstraint->second).linear[j].first, tmp);
593         }
594       }
595     }
596   }
597   inline void assemble(int entR, int typeR, const dataMat &value)
598   {
599     assemble(Dof(entR, typeR), value);
600   }
601   inline void assemble(MVertex *vR, int iCompR, int iFieldR,
602                        const dataMat &value)
603   {
604     assemble(vR->getNum(), Dof::createTypeWithTwoInts(iCompR, iFieldR), value);
605   }
606   virtual int sizeOfR() const
607   {
608     return _isParallel ? _localSize : unknown.size();
609   }
610   virtual int sizeOfF() const { return fixed.size(); }
611   virtual void systemSolve() { _current->systemSolve(); }
612   virtual void systemClear()
613   {
614     _current->zeroMatrix();
615     _current->zeroRightHandSide();
616   }
617   virtual inline void setCurrentMatrix(std::string name)
618   {
619     typename std::map<const std::string, linearSystem<dataMat> *>::iterator it =
620       _linearSystems.find(name);
621     if(it != _linearSystems.end())
622       _current = it->second;
623     else {
624       Msg::Error("Current matrix %s not found ", name.c_str());
625       throw;
626     }
627   }
628   virtual linearSystem<dataMat> *getLinearSystem(std::string &name)
629   {
630     typename std::map<const std::string, linearSystem<dataMat> *>::iterator it =
631       _linearSystems.find(name);
632     if(it != _linearSystems.end())
633       return it->second;
634     else
635       return 0;
636   }
637   virtual inline void
638   setLinearConstraint(Dof key, DofAffineConstraint<dataVec> &affineconstraint)
639   {
640     constraints[key] = affineconstraint;
641     // constraints.insert(std::make_pair(key, affineconstraint));
642   }
643 
644   virtual inline bool
645   getLinearConstraint(Dof key, DofAffineConstraint<dataVec> &affineconstraint)
646   {
647     typename std::map<Dof, DofAffineConstraint<dataVec> >::const_iterator it =
648       constraints.find(key);
649     if(it != constraints.end()) {
650       affineconstraint = it->second;
651       return true;
652     }
653     return false;
654   }
655 
656   virtual inline void assembleLinConst(const Dof &R, const Dof &C,
657                                        const dataMat &value)
658   {
659     auto itR = unknown.find(R);
660     if(itR != unknown.end()) {
661       typename std::map<Dof, DofAffineConstraint<dataVec> >::iterator
662         itConstraint;
663       itConstraint = constraints.find(C);
664       if(itConstraint != constraints.end()) {
665         dataMat tmp(value);
666         for(unsigned i = 0; i < (itConstraint->second).linear.size(); i++) {
667           dofTraits<T>::gemm(tmp, (itConstraint->second).linear[i].second,
668                              value, 1, 0);
669           assemble(R, (itConstraint->second).linear[i].first, tmp);
670         }
671         dataMat tmp2(value);
672         dofTraits<T>::gemm(tmp2, value, itConstraint->second.shift, -1, 0);
673         _current->addToRightHandSide(itR->second, tmp2);
674       }
675     }
676     else { // test function ; (no shift ?)
677       typename std::map<Dof, DofAffineConstraint<dataVec> >::iterator
678         itConstraint;
679       itConstraint = constraints.find(R);
680       if(itConstraint != constraints.end()) {
681         dataMat tmp(value);
682         for(unsigned i = 0; i < (itConstraint->second).linear.size(); i++) {
683           dofTraits<T>::gemm(tmp, itConstraint->second.linear[i].second, value,
684                              1, 0);
685           assemble((itConstraint->second).linear[i].first, C, tmp);
686         }
687       }
688     }
689   }
690   virtual void getFixedDof(std::vector<Dof> &R)
691   {
692     R.clear();
693     R.reserve(fixed.size());
694     typename std::map<Dof, dataVec>::iterator it;
695     for(it = fixed.begin(); it != fixed.end(); ++it) {
696       R.push_back(it->first);
697     }
698   }
699   virtual void getFixedDof(std::set<Dof> &R)
700   {
701     R.clear();
702     typename std::map<Dof, dataVec>::iterator it;
703     for(it = fixed.begin(); it != fixed.end(); ++it) {
704       R.insert(it->first);
705     }
706   }
707 
708   virtual int getDofNumber(const Dof &ky)
709   {
710     Dof key = ky;
711     {
712       auto it = associatedWith.find(ky);
713       if (it != associatedWith.end())key = it->second;
714     }
715 
716     auto it = unknown.find(key);
717     if(it == unknown.end()) {
718       return -1;
719     }
720     else
721       return it->second;
722   }
723 
724   virtual void clearAllLineConstraints() { constraints.clear(); }
725 
726   std::map<Dof, DofAffineConstraint<dataVec> > &getAllLinearConstraints()
727   {
728     return constraints;
729   };
730 };
731 
732 #endif
733