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
#ifndef GLOBAL_GEODESIC_FE_ASSEMBLER_HH
#define GLOBAL_GEODESIC_FE_ASSEMBLER_HH
#include <dune/istl/bcrsmatrix.hh>
#include <dune/common/fmatrix.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/matrix.hh>
#include "localgeodesicfestiffness.hh"
/** \brief A global FE assembler for problems involving functions that map into non-Euclidean spaces
*/
template <class GridView, class TargetSpace>
class GeodesicFEAssembler {
typedef typename GridView::template Codim<0>::Entity EntityType;
typedef typename GridView::template Codim<0>::EntityPointer EntityPointer;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
//! Dimension of the grid.
enum { gridDim = GridView::dimension };
//! Dimension of a tangent space
enum { blocksize = TargetSpace::TangentVector::size };
//!
typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
const GridView gridView_;
const LocalGeodesicFEStiffness<GridView,TargetSpace>* localStiffness_;
public:
/** \brief Constructor for a given grid */
GeodesicFEAssembler(const GridView& gridView) :
gridView_(gridView)
{}
/** \brief Assemble the tangent stiffness matrix
*/
void assembleMatrix(const std::vector<TargetSpace>& sol,
Dune::BCRSMatrix<MatrixBlock>& matrix) const;
/** \brief Assemble the gradient */
void assembleGradient(const std::vector<TargetSpace>& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const;
/** \brief Assemble the gradient using a finite difference approximation */
void assembleGradientFD(const std::vector<TargetSpace>& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const;
/** \brief Compute the energy of a deformation state */
double computeEnergy(const std::vector<TargetSpace>& sol) const;
protected:
void getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const;
}; // end class
template <class GridView, class TargetSpace>
void GeodesicFEAssembler<GridView,TargetSpace>::
getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
{
const typename GridView::IndexSet& indexSet = gridView_->indexSet();
int n = gridView_->size(gridDim);
nb.resize(n, n);
ElementIterator it = gridView_->template begin<0>();
ElementIterator endit = gridView_->template end<0> ();
for (; it!=endit; ++it) {
for (int i=0; i<it->template count<gridDim>(); i++) {
for (int j=0; j<it->template count<gridDim>(); j++) {
int iIdx = indexSet.subIndex(*it,i,gridDim);
int jIdx = indexSet.subIndex(*it,j,gridDim);
nb.add(iIdx, jIdx);
}
}
}
}
template <class GridView, class TargetSpace>
double GeodesicFEAssembler<GridView, TargetSpace>::
computeEnergy(const std::vector<TargetSpace>& sol) const
{
double energy = 0;
const typename GridView::IndexSet& indexSet = gridView_.indexSet();
if (sol.size()!=indexSet.size(gridDim))
DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
std::vector<TargetSpace> localSolution;
ElementIterator it = gridView_.template begin<0>();
ElementIterator endIt = gridView_.template end<0>();
// Loop over all elements
for (; it!=endIt; ++it) {
localSolution.resize(it->template count<gridDim>());
for (int i=0; i<it->template count<gridDim>(); i++)
localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
energy += localStiffness_->energy(*it, localSolution);
}
return energy;
}
#endif