From eca64c9339854d54c98433c7ee6c34dba5c77c90 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Wed, 14 Jan 2015 20:33:37 +0000
Subject: [PATCH] Read the initial iterate from a Python string, instead of
hardcoding it
This is much more flexible and *much* easier to use.
Along the way, this patch scraps several old hard-coded Dirichlet boundary
values. I don't even remember what they were good for. If they are needed
again they should also be reimplemented in Python.
[[Imported from SVN: r10007]]
---
src/harmonicmaps.cc | 73 ++++++---------------------------------------
1 file changed, 9 insertions(+), 64 deletions(-)
diff --git a/src/harmonicmaps.cc b/src/harmonicmaps.cc
index 9f089b8b..9d0697b8 100644
--- a/src/harmonicmaps.cc
+++ b/src/harmonicmaps.cc
@@ -25,6 +25,7 @@
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functions/vtkbasisgridfunction.hh>
+#include <dune/fufem/functiontools/basisinterpolator.hh>
#include <dune/fufem/dunepython.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
@@ -164,71 +165,17 @@ int main (int argc, char *argv[]) try
// Initial iterate
// //////////////////////////
- FieldVector<double,3> yAxis(0);
- yAxis[1] = 1;
-
- GridType::LeafGridView::Codim<dim>::Iterator vIt = grid.leafbegin<dim>();
- GridType::LeafGridView::Codim<dim>::Iterator vEndIt = grid.leafend<dim>();
-
- for (; vIt!=vEndIt; ++vIt) {
- int idx = grid.leafIndexSet().index(*vIt);
+ typedef P1NodalBasis<typename GridType::LeafGridView,double> FEBasis;
+ FEBasis feBasis(grid.leafGridView());
- TargetSpace::CoordinateType v;
+ std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialIterate") + std::string(")");
+ PythonFunction<FieldVector<double,dim>, TargetSpace::CoordinateType > pythonInitialIterate(Python::evaluate(lambda));
- FieldVector<double,dim> pos = vIt->geometry().corner(0);
- FieldVector<double,3> axis;
- axis[0] = pos[0]; axis[1] = pos[1]; axis[2] = 1;
- Rotation<double,3> rotation(axis, pos.two_norm()*M_PI*1.5);
+ std::vector<TargetSpace::CoordinateType> v;
+ Functions::interpolate(feBasis, v, pythonInitialIterate);
- //dirichletNodes[idx] = pos[0]<0.01 || pos[0] > 0.99;
-
- if (dirichletNodes[idx][0]) {
-// FieldMatrix<double,3,3> rMat;
-// rotation.matrix(rMat);
-// v = rMat[2];
-#ifdef UNITVECTOR3
- v[0] = std::sin(pos[0]*M_PI);
- v[1] = 0;
- v[2] = std::cos(pos[0]*M_PI);
-#endif
-#ifdef UNITVECTOR3
- v[0] = 1;
- v[1] = 0;
- v[2] = 0;
-#endif
-#if defined UNITVECTOR4 || defined ROTATION3
- FieldMatrix<double,3,3> rMat;
- rotation.matrix(rMat);
- v[0] = rMat[2][0];
- v[1] = rMat[2][1];
- v[2] = rMat[2][2];
- v[3] = 0;
-#endif
-#ifdef UNITVECTOR2
- v[0] = std::sin(pos[0]*M_PI);
- v[1] = std::cos(pos[0]*M_PI);
-#endif
-#if defined ROTATION2 || defined REALTUPLE1
- v[0] = pos[0]/*M_PI*/;
-#endif
- } else {
-#if defined UNITVECTOR2 || defined UNITVECTOR3 || defined UNITVECTOR4 || defined ROTATION3
- v = 0;
- v[0] = 1;
-#endif
-#if defined ROTATION2 || defined REALTUPLE1
- v[0] = 0.5*M_PI;
-#endif
- }
-
-#if defined UNITVECTOR2 || defined UNITVECTOR3 || defined UNITVECTOR4 || defined ROTATION3 || defined REALTUPLE1
- x[idx] = TargetSpace(v);
-#endif
-#if defined ROTATION2
- x[idx] = v[0];
-#endif
-
- }
+ for (size_t i=0; i<x.size(); i++)
+ x[i] = v[i];
// backup for error measurement later
SolutionType initialIterate = x;
@@ -236,8 +183,6 @@ int main (int argc, char *argv[]) try
// ////////////////////////////////////////////////////////////
// Create an assembler for the Harmonic Energy Functional
// ////////////////////////////////////////////////////////////
- typedef P1NodalBasis<typename GridType::LeafGridView,double> FEBasis;
- FEBasis feBasis(grid.leafGridView());
GFE::ChiralSkyrmionEnergy<GridType::LeafGridView, FEBasis::LocalFiniteElement, double> chiralSkyrmionEnergy;
HarmonicEnergyLocalStiffness<GridType::LeafGridView, FEBasis::LocalFiniteElement, TargetSpace> harmonicEnergyLocalStiffness;
--
GitLab