1 #ifndef DUNE_PDELAB_MFDCOMMON_HH
2 #define DUNE_PDELAB_MFDCOMMON_HH
6 #include <dune/common/fvector.hh>
7 #include <dune/geometry/referenceelements.hh>
14 template<
class ctype,
int dim>
17 typedef FieldVector<ctype,dim>
Vector;
30 template<
class Entity>
31 void init(
const Entity& ent)
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);
48 template<
class Intersection>
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();
64 for (
unsigned int i = 0; i <
dim; ++i)
65 R.push_back((face_center[i]-
center[i]) * face_area);
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]);
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);
91 template<
class ctype,
int dim>
99 unsigned int R_max = cell.
R.size();
104 for(
unsigned int i = 0; i <
dim; ++i)
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];
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)
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];
138 mutable std::vector<ctype> R_tilde;
139 mutable std::vector<ctype> sp;
142 template<
class ctype,
int dim>
151 const FieldMatrix<ctype,dim,dim>& K, std::vector<ctype>& W)
const
155 W.resize(this->
W_indep.size());
157 for (
int i = 0; i <
dim; ++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)
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];
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