From bf9d46ec58b6dab157412e6213bb110bcb34d14d Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Mon, 22 Aug 2011 10:27:52 +0000
Subject: [PATCH] Use Functions::interpolate to sample the analytically given
Dirichlet values
[[Imported from SVN: r7620]]
---
harmonicmaps-eoc.cc | 69 ++++++++++++++++++++++-----------------------
1 file changed, 33 insertions(+), 36 deletions(-)
diff --git a/harmonicmaps-eoc.cc b/harmonicmaps-eoc.cc
index 74ff8928..5d6c4b0f 100644
--- a/harmonicmaps-eoc.cc
+++ b/harmonicmaps-eoc.cc
@@ -18,6 +18,7 @@
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
+#include <dune/fufem/functiontools/basisinterpolator.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/norms/h1seminorm.hh>
@@ -39,6 +40,23 @@ const int blocksize = TargetSpace::TangentVector::dimension;
using namespace Dune;
using std::string;
+struct DirichletFunction
+ : public Dune::VirtualFunction<FieldVector<double,dim>, FieldVector<double,3> >
+{
+ void evaluate(const FieldVector<double, dim>& x, FieldVector<double,3>& out) const {
+
+ FieldVector<double,3> axis;
+ axis[0] = x[0]; axis[1] = x[1]; axis[2] = 1;
+ Rotation<3,double> rotation(axis, x.two_norm()*M_PI*3);
+
+ FieldMatrix<double,3,3> rMat;
+ rotation.matrix(rMat);
+ out = rMat[2];
+
+ }
+};
+
+
template <class GridType>
void solve (const shared_ptr<GridType>& grid,
SolutionType& x,
@@ -68,44 +86,23 @@ void solve (const shared_ptr<GridType>& grid,
// Initial solution
// //////////////////////////
- x.resize(grid->size(dim));
-
- FieldVector<double,3> yAxis(0);
- yAxis[1] = 1;
-
- typename GridType::LeafGridView::template Codim<dim>::Iterator vIt = grid->template leafbegin<dim>();
- typename GridType::LeafGridView::template Codim<dim>::Iterator vEndIt = grid->template leafend<dim>();
-
- for (; vIt!=vEndIt; ++vIt) {
- int idx = grid->leafIndexSet().index(*vIt);
-
- FieldVector<double,3> v;
- FieldVector<double,dim> pos = vIt->geometry().corner(0);
- FieldVector<double,3> axis;
- axis[0] = pos[0]; axis[1] = pos[1]; axis[2] = 1;
- Rotation<3,double> rotation(axis, pos.two_norm()*M_PI*3);
-
- if (dirichletNodes[idx][0]) {
-#if 1
- FieldMatrix<double,3,3> rMat;
- rotation.matrix(rMat);
- v = rMat[2];
-#else
- v[0] = std::sin(pos[0]*M_PI);
- v[1] = 0;
- v[2] = std::cos(pos[0]*M_PI);
-#endif
- } else {
- v[0] = 1;
- v[1] = 0;
- v[2] = 0;
- }
-
- x[idx] = v;
-
- }
+ typedef P1NodalBasis<typename GridType::LeafGridView,double> P1Basis;
+ P1Basis p1Basis(grid->leafView());
+
+ x.resize(p1Basis.size());
+ BlockVector<FieldVector<double,3> > dirichletFunctionValues;
+ DirichletFunction dirichletFunction;
+ Functions::interpolate(p1Basis, dirichletFunctionValues, dirichletFunction);
+ FieldVector<double,3> innerValue;
+ innerValue[0] = 1;
+ innerValue[1] = 0;
+ innerValue[2] = 0;
+
+ for (size_t i=0; i<x.size(); i++)
+ x[i] = (dirichletNodes[i][0]) ? dirichletFunctionValues[i] : innerValue;
+
// ////////////////////////////////////////////////////////////
// Create an assembler for the Harmonic Energy Functional
// ////////////////////////////////////////////////////////////
--
GitLab