dune-pdelab  2.0.0
mfdcommon.hh
Go to the documentation of this file.
1 #ifndef DUNE_PDELAB_MFDCOMMON_HH
2 #define DUNE_PDELAB_MFDCOMMON_HH
3 
4 #include <iomanip>
5 #include <vector>
6 #include <dune/common/fvector.hh>
7 #include <dune/geometry/referenceelements.hh>
8 
9 namespace Dune
10 {
11  namespace PDELab
12  {
13  // Cache for cell information used by mimetic finite difference mthods
14  template<class ctype,int dim>
16  {
17  typedef FieldVector<ctype,dim> Vector;
18 
20  ctype volume;
21  unsigned int num_faces;
22 
23  std::vector<ctype> face_areas;
24  std::vector<Vector> face_centers;
25  std::vector<bool> boundary;
26  std::vector<ctype> R;
27  std::vector<ctype> N;
28 
29  // initialise MimeticCellProperties for given entity
30  template<class Entity>
31  void init(const Entity& ent)
32  {
33  const typename Entity::Geometry& cgeo = ent.geometry();
34  const FieldVector<ctype,dim> center_local
35  = ReferenceElements<ctype,dim>::general(cgeo.type()).position(0,0);
36  center = cgeo.global(center_local);
37  volume = cgeo.volume();
38  num_faces = 0;
39 
40  face_areas.clear();
41  face_centers.clear();
42  boundary.clear();
43  R.clear();
44  N.clear();
45  }
46 
47  // Append information on face described by the intersection object
48  template<class Intersection>
49  void add_face(const Intersection& is)
50  {
51  const typename Intersection::Geometry& fgeo = is.geometry();
52  const FieldVector<ctype,dim-1> face_center_local
53  = ReferenceElements<ctype,dim-1>::general(fgeo.type()).position(0,0);
54  const Vector face_center = fgeo.global(face_center_local);
55  const ctype face_area = fgeo.volume();
56 
57  face_areas.push_back(face_area);
58  face_centers.push_back(face_center);
59  boundary.push_back(is.boundary());
60 
61  ++num_faces;
62 
63  // construct matrix R
64  for (unsigned int i = 0; i < dim; ++i)
65  R.push_back((face_center[i]-center[i]) * face_area);
66 
67  // construct matrix N
68  Vector face_normal = is.unitOuterNormal(face_center_local);
69  for (unsigned int i = 0; i < dim; ++i)
70  N.push_back(face_normal[i]);
71  }
72 
73  void print() const
74  {
75  std::cout.precision(3);
76  std::cout << "Cell center (" << center[0];
77  for (unsigned int i = 1; i < dim; ++i)
78  std::cout << ", " << center[i];
79  std::cout << "), volume " << volume << std::endl;
80  std::cout << "face areas " << std::setw(6) << face_areas[0];
81  for (unsigned i = 1; i < face_areas.size(); ++i)
82  std::cout << ", " << std::setw(6) << face_areas[i];
83  std::cout << "\n is boundary " << std::setw(4) << boundary[0];
84  for (unsigned i = 1; i < boundary.size(); ++i)
85  std::cout << ", " << std::setw(4) << boundary[i];
86  std::cout << std::endl;
87  std::cout.precision(6);
88  }
89  };
90 
91  template<class ctype,int dim>
93  {
94  protected:
95  mutable std::vector<ctype> W_indep;
96 
98  {
99  unsigned int R_max = cell.R.size();
100 
101  // orthonormalize R
102  R_tilde = cell.R;
103  sp.reserve(dim - 1);
104  for(unsigned int i = 0; i < dim; ++i)
105  {
106  // orthogonalize
107  sp.clear(); sp.resize(i, 0.0);
108  for(unsigned int j = 0; j < i; ++j)
109  for(unsigned int k = 0; k < R_max; k += dim)
110  sp[j] += R_tilde[k+i] * R_tilde[k+j];
111  for(unsigned int j = 0; j < i; ++j)
112  for(unsigned int k = 0; k < R_max; k += dim)
113  R_tilde[k+i] -= sp[j]*R_tilde[k+j];
114 
115  // normalize
116  ctype norm = 0.0;
117  for(unsigned int k = i; k < R_max; k += dim)
118  norm += R_tilde[k] * R_tilde[k];
119  norm = 1.0 / std::sqrt(norm);
120  for(unsigned int k = i; k < R_max; k += dim)
121  R_tilde[k] *= norm;
122  }
123 
124  // construct the part of matrix W that is independent of K
125  W_indep.clear();
126  W_indep.resize(cell.num_faces * cell.num_faces, 0.0);
127  for (unsigned int i = 0; i < W_indep.size(); i += cell.num_faces+1)
128  W_indep[i] = 1.0;
129  for (unsigned int i = 0, m = 0; i < R_max; i += dim)
130  for (unsigned int j = 0; j < R_max; j += dim, ++m)
131  for (unsigned int k = 0; k < dim; ++k)
132  W_indep[m] -= R_tilde[i+k] * R_tilde[j+k];
133  }
134 
135  private:
136  // local variables of build_W_indep; defined here only for
137  // performance reasons (no need to reallocate each time)
138  mutable std::vector<ctype> R_tilde;
139  mutable std::vector<ctype> sp;
140  };
141 
142  template<class ctype,int dim>
143  class MimeticBrezziW : public MimeticBrezziWBase<ctype,dim>
144  {
145  const ctype alpha;
146 
147  public:
148  explicit MimeticBrezziW(ctype alpha_ = 1.0) : alpha(alpha_) {}
149 
151  const FieldMatrix<ctype,dim,dim>& K, std::vector<ctype>& W) const
152  {
153  this->build_W_indep(cell);
154 
155  W.resize(this->W_indep.size());
156  ctype u = 0.0;
157  for (int i = 0; i < dim; ++i)
158  u += K[i][i];
159  u *= 2.0 * alpha / dim;
160  ctype inv_volume = 1.0 / cell.volume;
161  for (unsigned int i = 0, m = 0; i < cell.num_faces; ++i)
162  for (unsigned int j = 0; j < cell.num_faces; ++j, ++m)
163  {
164  W[m] = u * this->W_indep[m];
165  for (unsigned int k = 0; k < dim; ++k)
166  for (unsigned int l = 0; l < dim; ++l)
167  W[m] += cell.N[i*dim+k] * K[k][l] * cell.N[j*dim+l];
168  W[m] *= cell.face_areas[i] * cell.face_areas[j] * inv_volume;
169  }
170  }
171  };
172  }
173 }
174 
175 #endif
std::vector< ctype > face_areas
Definition: mfdcommon.hh:23
Definition: mfdcommon.hh:143
void add_face(const Intersection &is)
Definition: mfdcommon.hh:49
ctype volume
Definition: mfdcommon.hh:20
std::vector< ctype > R
Definition: mfdcommon.hh:26
static const int dim
Definition: adaptivity.hh:82
std::vector< Vector > face_centers
Definition: mfdcommon.hh:24
Vector center
Definition: mfdcommon.hh:19
unsigned int num_faces
Definition: mfdcommon.hh:21
Definition: mfdcommon.hh:15
std::vector< ctype > W_indep
Definition: mfdcommon.hh:95
MimeticBrezziW(ctype alpha_=1.0)
Definition: mfdcommon.hh:148
void build_W_indep(const MimeticCellProperties< ctype, dim > &cell) const
Definition: mfdcommon.hh:97
void init(const Entity &ent)
Definition: mfdcommon.hh:31
void build_W(const MimeticCellProperties< ctype, dim > &cell, const FieldMatrix< ctype, dim, dim > &K, std::vector< ctype > &W) const
Definition: mfdcommon.hh:150
Definition: mfdcommon.hh:92
std::vector< ctype > N
Definition: mfdcommon.hh:27
std::vector< bool > boundary
Definition: mfdcommon.hh:25
FieldVector< ctype, dim > Vector
Definition: mfdcommon.hh:17
void print() const
Definition: mfdcommon.hh:73