From 0b884ab37d110c8cd672fabc20316676e8753fa3 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 18 Mar 2010 18:19:55 +0000
Subject: [PATCH] globally refine a grid and transfer a given function to the
next level. Works for all TargetSpaces and replace rodRefine, which did the
same only for RigidBodyMotions
[[Imported from SVN: r5771]]
---
src/Makefile.am | 2 +-
src/geodesicfefunctionadaptor.hh | 85 ++++++++++++++++++++++++++++++++
src/rodrefine.hh | 78 -----------------------------
3 files changed, 86 insertions(+), 79 deletions(-)
create mode 100644 src/geodesicfefunctionadaptor.hh
delete mode 100644 src/rodrefine.hh
diff --git a/src/Makefile.am b/src/Makefile.am
index 3a455585..36a2c8cd 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -6,7 +6,7 @@ AM_CPPFLAGS = @AM_CPPFLAGS@ -I$(top_srcdir)/..
srcincludedir = $(includedir)/dune/common
srcinclude_HEADERS = averagedistanceassembler.hh localgeodesicfefunction.hh \
- quaternion.hh rodrefine.hh averageinterface.hh localstiffness.hh \
+ quaternion.hh geodesicfefunctionadaptor.hh averageinterface.hh localstiffness.hh \
localgeodesicfestiffness.hh riemanniantrsolver.hh rodwriter.hh \
geodesicdifference.hh makestraightrod.hh rigidbodymotion.hh \
rotation.hh geodesicfeassembler.hh maxnormtrustregion.hh \
diff --git a/src/geodesicfefunctionadaptor.hh b/src/geodesicfefunctionadaptor.hh
new file mode 100644
index 00000000..b9340a67
--- /dev/null
+++ b/src/geodesicfefunctionadaptor.hh
@@ -0,0 +1,85 @@
+#ifndef GEODESIC_FE_FUNCTION_ADAPTOR_HH
+#define GEODESIC_FE_FUNCTION_ADAPTOR_HH
+
+#include <vector>
+#include <map>
+
+#include "localgeodesicfefunction.hh"
+
+template <class GridType, class TargetSpace>
+void geodesicFEFunctionAdaptor(GridType& grid, std::vector<TargetSpace>& x)
+{
+ const int dim = GridType::dimension;
+
+ typedef typename GridType::template Codim<0>::LeafIterator ElementIterator;
+ typedef typename GridType::template Codim<dim>::LeafIterator VertexIterator;
+
+ // /////////////////////////////////////////////////////
+ // Save leaf p1 data in a map
+ // /////////////////////////////////////////////////////
+
+ const typename GridType::Traits::LocalIdSet& idSet = grid.localIdSet();
+ const typename GridType::Traits::LeafIndexSet& indexSet = grid.leafIndexSet();
+
+ std::map<typename GridType::Traits::LocalIdSet::IdType, TargetSpace> dofMap;
+
+ VertexIterator vIt = grid.template leafbegin<dim>();
+ VertexIterator vEndIt = grid.template leafend<dim>();
+
+ for (; vIt!=vEndIt; ++vIt)
+ dofMap.insert(std::make_pair(idSet.id(*vIt), x[indexSet.index(*vIt)]));
+
+
+
+ // /////////////////////////////////////////////////////
+ // Globally refine the grid
+ // /////////////////////////////////////////////////////
+
+ grid.globalRefine(1);
+
+
+ // /////////////////////////////////////////////////////
+ // Restore and interpolate the data
+ // /////////////////////////////////////////////////////
+
+ x.resize(grid.size(dim));
+
+ ElementIterator eIt = grid.template leafbegin<0>();
+ ElementIterator eEndIt = grid.template leafend<0>();
+
+ for (; eIt!=eEndIt; ++eIt) {
+
+ // Set up a local gfe function on the father element
+ std::vector<TargetSpace> coefficients(eIt->father()->template count<dim>());
+
+ for (int i=0; i<eIt->father()->template count<dim>(); i++)
+ coefficients[i] = dofMap.find(idSet.subId(*eIt->father(),i,dim))->second;
+
+ LocalGeodesicFEFunction<dim,double,TargetSpace> fatherFunction(coefficients);
+
+ // The embedding of this element into the father geometry
+ const typename GridType::template Codim<0>::LocalGeometry& geometryInFather = eIt->geometryInFather();
+
+ for (int i=0; i<eIt->template count<dim>(); i++) {
+
+ if (dofMap.find(idSet.subId(*eIt,i,dim)) != dofMap.end()) {
+
+ // If the vertex exists on the coarser level we take the value from there.
+ // That should be faster and more accurate than interpolating
+ x[indexSet.subIndex(*eIt,i,dim)] = dofMap[idSet.subId(*eIt,i,dim)];
+
+ } else {
+
+ // Interpolate
+ x[indexSet.subIndex(*eIt,i,dim)] = fatherFunction.evaluate(geometryInFather.corner(i));
+
+ }
+
+ }
+
+ }
+
+
+}
+
+#endif
diff --git a/src/rodrefine.hh b/src/rodrefine.hh
deleted file mode 100644
index d847035c..00000000
--- a/src/rodrefine.hh
+++ /dev/null
@@ -1,78 +0,0 @@
-#ifndef ROD_REFINE_HH
-#define ROD_REFINE_HH
-
-#include <vector>
-#include <map>
-
-#include "rigidbodymotion.hh"
-
-template <class GridType>
-void globalRodRefine(GridType& grid, std::vector<RigidBodyMotion<3> >& x)
-{
- // Should be a static assertion
- assert(GridType::dimension == 1);
-
- typedef typename GridType::template Codim<0>::LeafIterator ElementIterator;
- typedef typename GridType::template Codim<1>::LeafIterator VertexIterator;
-
- // /////////////////////////////////////////////////////
- // Save leaf p1 data in a map
- // /////////////////////////////////////////////////////
-
- const typename GridType::Traits::LocalIdSet& idSet = grid.localIdSet();
- const typename GridType::Traits::LeafIndexSet& indexSet = grid.leafIndexSet();
-
- std::map<typename GridType::Traits::LocalIdSet::IdType, RigidBodyMotion<3> > dofMap;
-
- VertexIterator vIt = grid.template leafbegin<1>();
- VertexIterator vEndIt = grid.template leafend<1>();
-
- for (; vIt!=vEndIt; ++vIt)
- dofMap.insert(std::make_pair(idSet.id(*vIt), x[indexSet.index(*vIt)]));
-
-
-
- // /////////////////////////////////////////////////////
- // Globally refine the grid
- // /////////////////////////////////////////////////////
-
- grid.globalRefine(1);
-
-
- // /////////////////////////////////////////////////////
- // Restore and interpolate the data
- // /////////////////////////////////////////////////////
-
- x.resize(grid.size(1));
-
- ElementIterator eIt = grid.template leafbegin<0>();
- ElementIterator eEndIt = grid.template leafend<0>();
-
- for (; eIt!=eEndIt; ++eIt) {
-
- for (int i=0; i<2; i++) {
-
- if (dofMap.find(idSet.subId(*eIt,i,1)) != dofMap.end()) {
-
- x[indexSet.subIndex(*eIt,i,1)] = dofMap[idSet.subId(*eIt,i,1)];
-
- } else {
- // Interpolate
- RigidBodyMotion<3> p0 = dofMap[idSet.subId(*eIt->father(),0,1)];
- RigidBodyMotion<3> p1 = dofMap[idSet.subId(*eIt->father(),1,1)];
-
- x[indexSet.subIndex(*eIt,i,1)].r = (p0.r + p1.r);
- x[indexSet.subIndex(*eIt,i,1)].r *= 0.5;
- x[indexSet.subIndex(*eIt,i,1)].q
- = Rotation<3,double>::interpolate(p0.q, p1.q, 0.5);
-
- }
-
- }
-
- }
-
-
-}
-
-#endif
--
GitLab