From 7618fbe0e1357d146978e99a2f14a07101d72619 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Fri, 17 Apr 2009 15:02:37 +0000
Subject: [PATCH] add methods to assemble the gradient and the Hesse matrix
[[Imported from SVN: r4030]]
---
src/geodesicfeassembler.hh | 95 ++++++++++++++++++++++++++++++++++++--
1 file changed, 90 insertions(+), 5 deletions(-)
diff --git a/src/geodesicfeassembler.hh b/src/geodesicfeassembler.hh
index 542ca535..e860cded 100644
--- a/src/geodesicfeassembler.hh
+++ b/src/geodesicfeassembler.hh
@@ -29,7 +29,7 @@ class GeodesicFEAssembler {
const GridView gridView_;
- const LocalGeodesicFEStiffness<GridView,TargetSpace>* localStiffness_;
+ LocalGeodesicFEStiffness<GridView,TargetSpace>* localStiffness_;
public:
@@ -65,14 +65,14 @@ template <class GridView, class TargetSpace>
void GeodesicFEAssembler<GridView,TargetSpace>::
getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
{
- const typename GridView::IndexSet& indexSet = gridView_->indexSet();
+ const typename GridView::IndexSet& indexSet = gridView_.indexSet();
- int n = gridView_->size(gridDim);
+ int n = gridView_.size(gridDim);
nb.resize(n, n);
- ElementIterator it = gridView_->template begin<0>();
- ElementIterator endit = gridView_->template end<0> ();
+ ElementIterator it = gridView_.template begin<0>();
+ ElementIterator endit = gridView_.template end<0> ();
for (; it!=endit; ++it) {
@@ -93,6 +93,91 @@ getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
}
+template <class GridView, class TargetSpace>
+void GeodesicFEAssembler<GridView,TargetSpace>::
+assembleMatrix(const std::vector<TargetSpace>& sol,
+ Dune::BCRSMatrix<MatrixBlock>& matrix) const
+{
+ const typename GridView::IndexSet& indexSet = gridView_.indexSet();
+
+ Dune::MatrixIndexSet neighborsPerVertex;
+ getNeighborsPerVertex(neighborsPerVertex);
+
+ matrix = 0;
+
+ ElementIterator it = gridView_.template begin<0>();
+ ElementIterator endit = gridView_.template end<0> ();
+
+ for( ; it != endit; ++it ) {
+
+ const int numOfBaseFct = it->template count<gridDim>();
+
+ // Extract local solution
+ std::vector<TargetSpace> localSolution(numOfBaseFct);
+
+ for (int i=0; i<numOfBaseFct; i++)
+ localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
+
+ // setup matrix
+ localStiffness_->assemble(*it, localSolution);
+
+ // Add element matrix to global stiffness matrix
+ for(int i=0; i<numOfBaseFct; i++) {
+
+ int row = indexSet.subIndex(*it,i,gridDim);
+
+ for (int j=0; j<numOfBaseFct; j++ ) {
+
+ int col = indexSet.subIndex(*it,j,gridDim);
+ matrix[row][col] += localStiffness_->mat(i,j);
+
+ }
+ }
+
+ }
+
+}
+
+template <class GridView, class TargetSpace>
+void GeodesicFEAssembler<GridView,TargetSpace>::
+assembleGradient(const std::vector<TargetSpace>& sol,
+ Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
+{
+ const typename GridView::IndexSet& indexSet = gridView_.indexSet();
+
+ if (sol.size()!=gridView_.size(gridDim))
+ DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
+
+ grad.resize(sol.size());
+ grad = 0;
+
+ ElementIterator it = gridView_->template begin<0>();
+ ElementIterator endIt = gridView_->template end<0>();
+
+ // Loop over all elements
+ for (; it!=endIt; ++it) {
+
+ // A 1d grid has two vertices
+ const int nDofs = it->template count<gridDim>();
+
+ // Extract local solution
+ std::vector<TargetSpace> localSolution(nDofs);
+
+ for (int i=0; i<nDofs; i++)
+ localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
+
+ // Assemble local gradient
+ std::vector<Dune::FieldVector<double,blocksize> > localGradient(nDofs);
+
+ localStiffness_->assembleGradient(*it, localSolution, localGradient);
+
+ // Add to global gradient
+ for (int i=0; i<nDofs; i++)
+ grad[indexSet.subIndex(*it,i,gridDim)] += localGradient[i];
+
+ }
+
+}
template <class GridView, class TargetSpace>
--
GitLab