Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
#ifndef GFE_DIFFERENCE_NORM_SQUARED_HH
#define GFE_DIFFERENCE_NORM_SQUARED_HH
/** \file
\brief Helper method which computes the energy norm (squared) of the difference of
the fe functions on two different refinements of the same base grid.
*/
#include <dune/geometry/type.hh>
#include <dune/common/fvector.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/fufem/assemblers/localoperatorassembler.hh>
#include <dune/gfe/globalgeodesicfefunction.hh>
template <class Basis, class TargetSpace>
class GFEDifferenceNormSquared {
typedef typename Basis::GridView GridView;
typedef typename Basis::GridView::Grid GridType;
// Grid dimension
const static int dim = GridType::dimension;
/** \brief Check whether given local coordinates are contained in the reference element
\todo This method exists in the Dune grid interface! But we need the eps.
*/
static bool checkInside(const Dune::GeometryType& type,
const Dune::FieldVector<double, dim> &loc,
double eps)
{
switch (type.dim()) {
case 0: // vertex
return false;
case 1: // line
return -eps <= loc[0] && loc[0] <= 1+eps;
case 2:
if (type.isSimplex()) {
return -eps <= loc[0] && -eps <= loc[1] && (loc[0]+loc[1])<=1+eps;
} else if (type.isCube()) {
return -eps <= loc[0] && loc[0] <= 1+eps
&& -eps <= loc[1] && loc[1] <= 1+eps;
} else {
DUNE_THROW(Dune::GridError, "checkInside(): ERROR: Unknown type " << type << " found!");
}
case 3:
if (type.isSimplex()) {
return -eps <= loc[0] && -eps <= loc[1] && -eps <= loc[2]
&& (loc[0]+loc[1]+loc[2]) <= 1+eps;
} else if (type.isPyramid()) {
return -eps <= loc[0] && -eps <= loc[1] && -eps <= loc[2]
&& (loc[0]+loc[2]) <= 1+eps
&& (loc[1]+loc[2]) <= 1+eps;
} else if (type.isPrism()) {
return -eps <= loc[0] && -eps <= loc[1]
&& (loc[0]+loc[1])<= 1+eps
&& -eps <= loc[2] && loc[2] <= 1+eps;
} else if (type.isCube()) {
return -eps <= loc[0] && loc[0] <= 1+eps
&& -eps <= loc[1] && loc[1] <= 1+eps
&& -eps <= loc[2] && loc[2] <= 1+eps;
}else {
DUNE_THROW(Dune::GridError, "checkInside(): ERROR: Unknown type " << type << " found!");
}
default:
DUNE_THROW(Dune::GridError, "checkInside(): ERROR: Unknown type " << type << " found!");
}
}
/** \brief Coordinate function in one variable, constant in the others
This is used to extract the positions of the Lagrange nodes.
*/
struct CoordinateFunction
: public Dune::VirtualFunction<Dune::FieldVector<double,dim>, Dune::FieldVector<double,1> >
{
CoordinateFunction(int d)
: d_(d)
{}
void evaluate(const Dune::FieldVector<double, dim>& x, Dune::FieldVector<double,1>& out) const {
out[0] = x[d_];
}
//
int d_;
};
public:
static void interpolate(const GridType& sourceGrid,
const std::vector<TargetSpace>& source,
const GridType& targetGrid,
std::vector<TargetSpace>& target)
{
// Create a leaf function, which we need to be able to call 'evalall()'
Basis sourceBasis(sourceGrid.leafView());
GlobalGeodesicFEFunction<Basis,TargetSpace> sourceFunction(sourceBasis, source);
Basis targetBasis(targetGrid.leafView());
// ///////////////////////////////////////////////////////////////////////////////////////////
// Prolong the adaptive solution onto the uniform grid in order to make it comparable
// ///////////////////////////////////////////////////////////////////////////////////////////
target.resize(targetBasis.size());
// handle each dof only once
Dune::BitSetVector<1> handled(targetBasis.size(), false);
typename GridView::template Codim<0>::Iterator eIt = targetGrid.template leafbegin<0>();
typename GridView::template Codim<0>::Iterator eEndIt = targetGrid.template leafend<0>();
// Loop over all dofs by looping over all elements
for (; eIt!=eEndIt; ++eIt) {
const typename Basis::LocalFiniteElement& targetLFE = targetBasis.getLocalFiniteElement(*eIt);
// Generate position of the Lagrange nodes
std::vector<Dune::FieldVector<double,dim> > lagrangeNodes(targetLFE.localBasis().size());
for (int i=0; i<dim; i++) {
CoordinateFunction lFunction(i);
std::vector<Dune::FieldVector<double,1> > coordinates;
targetLFE.localInterpolation().interpolate(lFunction, coordinates);
for (size_t j=0; j<coordinates.size(); j++)
lagrangeNodes[j][i] = coordinates[j];
}
for (size_t i=0; i<targetLFE.localCoefficients().size(); i++) {
typename GridType::template Codim<0>::EntityPointer element(eIt);
Dune::FieldVector<double,dim> pos = lagrangeNodes[i];
if (handled[targetBasis.index(*eIt,i)][0])
continue;
handled[targetBasis.index(*eIt,i)] = true;
assert(checkInside(element->type(), pos, 1e-7));
// ////////////////////////////////////////////////////////////////////
// Get an element on the coarsest grid which contains the vertex and
// its local coordinates there
// ////////////////////////////////////////////////////////////////////
while (element->level() != 0){
pos = element->geometryInFather().global(pos);
element = element->father();
assert(checkInside(element->type(), pos, 1e-7));
}
// ////////////////////////////////////////////////////////////////////
// Find the corresponding coarse grid element on the adaptive grid.
// This is a linear algorithm, but we expect the coarse grid to be small.
// ////////////////////////////////////////////////////////////////////
Dune::LevelMultipleCodimMultipleGeomTypeMapper<GridType,Dune::MCMGElementLayout> uniformP0Mapper (targetGrid, 0);
Dune::LevelMultipleCodimMultipleGeomTypeMapper<GridType,Dune::MCMGElementLayout> adaptiveP0Mapper(sourceGrid, 0);
int coarseIndex = uniformP0Mapper.map(*element);
typename GridType::template Codim<0>::LevelIterator adaptEIt = sourceGrid.template lbegin<0>(0);
typename GridType::template Codim<0>::LevelIterator adaptEEndIt = sourceGrid.template lend<0>(0);
for (; adaptEIt!=adaptEEndIt; ++adaptEIt)
if (adaptiveP0Mapper.map(*adaptEIt) == coarseIndex)
break;
element = adaptEIt;
// ////////////////////////////////////////////////////////////////////////
// Find a corresponding point (not necessarily vertex) on the leaf level
// of the adaptive grid.
// ////////////////////////////////////////////////////////////////////////
while (!element->isLeaf()) {
typename GridType::template Codim<0>::Entity::HierarchicIterator hIt = element->hbegin(element->level()+1);
typename GridType::template Codim<0>::Entity::HierarchicIterator hEndIt = element->hend(element->level()+1);
Dune::FieldVector<double,dim> childPos;
assert(checkInside(element->type(), pos, 1e-7));
for (; hIt!=hEndIt; ++hIt) {
childPos = hIt->geometryInFather().local(pos);
if (checkInside(hIt->type(), childPos, 1e-7))
break;
}
assert(hIt!=hEndIt);
element = hIt;
pos = childPos;
}
// ////////////////////////////////////////////////////////////////////////
// Sample adaptive function
// ////////////////////////////////////////////////////////////////////////
sourceFunction.evaluateLocal(*element, pos,
target[targetBasis.index(*eIt,i)]);
}
}
}
template <int blocksize>
static double computeNormSquared(const GridType& grid,
const Dune::BlockVector<Dune::FieldVector<double,blocksize> >& x,
const LocalOperatorAssembler<GridType,
typename Basis::LocalFiniteElement,
typename Basis::LocalFiniteElement,
Dune::FieldMatrix<double,1,1> >* localStiffness)
{
// ///////////////////////////////////////////////////////////////////////////////////////////
// Compute the energy norm
// ///////////////////////////////////////////////////////////////////////////////////////////
double energyNormSquared = 0;
Basis basis(grid.leafView());
Dune::Matrix<Dune::FieldMatrix<double,1,1> > localMatrix;
typename GridType::template Codim<0>::LeafIterator eIt = grid.template leafbegin<0>();
typename GridType::template Codim<0>::LeafIterator eEndIt = grid.template leafend<0>();
for (; eIt!=eEndIt; ++eIt) {
size_t nLocalDofs = basis.getLocalFiniteElement(*eIt).localCoefficients().size();
localMatrix.setSize(nLocalDofs,nLocalDofs);
localStiffness->assemble(*eIt, localMatrix,
basis.getLocalFiniteElement(*eIt),
basis.getLocalFiniteElement(*eIt));
for (size_t i=0; i<nLocalDofs; i++)
for (size_t j=0; j<nLocalDofs; j++)
energyNormSquared += localMatrix[i][j][0][0] * (x[basis.index(*eIt,i)] * x[basis.index(*eIt,j)]);
}
return energyNormSquared;
}
static double compute(const GridType& uniformGrid,
const std::vector<TargetSpace>& uniformSolution,
const GridType& adaptiveGrid,
const std::vector<TargetSpace>& adaptiveSolution,
const LocalOperatorAssembler<GridType,
typename Basis::LocalFiniteElement,
typename Basis::LocalFiniteElement,
Dune::FieldMatrix<double,1,1> >* localStiffness)
{
std::vector<TargetSpace> uniformAdaptiveSolution;
interpolate(adaptiveGrid, adaptiveSolution, uniformGrid, uniformAdaptiveSolution);
// Compute difference in the embedding space
Dune::BlockVector<typename TargetSpace::CoordinateType> difference(uniformSolution.size());
for (size_t i=0; i<difference.size(); i++)
difference[i] = uniformAdaptiveSolution[i].globalCoordinates() - uniformSolution[i].globalCoordinates();
//std::cout << "new difference:\n" << difference << std::endl;
// std::cout << "new projected:\n" << std::endl;
// for (size_t i=0; i<difference.size(); i++)
// std::cout << uniformAdaptiveSolution[i].globalCoordinates() << std::endl;
return computeNormSquared(uniformGrid, difference, localStiffness);
}
};
#endif