added: functions to interpolate (or rather, for now; approximate) a Field for ASMs

git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@2701 e10b68d5-8a6e-419e-a041-bce267b0401d

src/ASM/ASMbase.h

@@ -28,6 +28,7 @@ typedef MPCSet::const_iterator MPCIter; //!< Iterator over an MPC equation set
 
 struct TimeDomain;
 class ElementBlock;
+class Field;
 class GlobalIntegral;
 class IntegrandBase;
 class Integrand;
@@ -413,6 +414,11 @@ public:
   virtual bool evaluate(const ASMbase* basis, const Vector& locVec,
                         Vector& vec) const { return false; }
 
+  //! \brief Evaluates and interpolates a field over a given geometry.
+  //! \param[in] field The field to evaluate
+  //! \param[out] vec The obtained coefficients after interpolation
+  virtual bool evaluate(const Field* field, Vector& vec) const { return false; }
+
   //! \brief Evaluates the secondary solution field at all visualization points.
   //! \param[out] sField Solution field
   //! \param[in] integrand Object with problem-specific data and methods

src/ASM/ASMs2D.h

@@ -341,6 +341,11 @@ public:
   virtual bool evaluate(const ASMbase* basis, const Vector& locVec,
                         Vector& vec) const;
 
+  //! \copydoc ASMbase::evaluate(const Field*,Vector&)
+  //! \details Note a VDSA is used as the regular interpolation method in
+  //!          GoTools only works with uniform knots.
+  virtual bool evaluate(const Field* field, Vector& vec) const;
+
   //! \brief Evaluates the secondary solution field at all visualization points.
   //! \param[out] sField Solution field
   //! \param[in] integrand Object with problem-specific data and methods

src/ASM/ASMs2Drecovery.C

@@ -15,6 +15,8 @@
 #include "GoTools/geometry/SurfaceInterpolator.h"
 
 #include "ASMs2D.h"
+#include "FiniteElement.h"
+#include "Field.h"
 #include "IntegrandBase.h"
 #include "CoordinateMapping.h"
 #include "GaussQuadrature.h"
@@ -456,6 +458,56 @@ Go::SplineSurface* ASMs2D::scRecovery (const IntegrandBase& integrand) const
 
 #include "ASMs2DInterpolate.C" // TODO: inline these methods instead...
 
+
+bool ASMs2D::evaluate (const Field* field, Vector& vec) const
+{
+  // Compute parameter values of the result sampling points (Greville points)
+  RealArray gpar[2];
+  for (int dir = 0; dir < 2; dir++)
+    if (!this->getGrevilleParameters(gpar[dir],dir))
+      return false;
+
+  // Evaluate the result field at all sampling points.
+  // Note: it is here assumed that *basis and *this have spline bases
+  // defined over the same parameter domain.
+  Vector sValues(gpar[0].size()*gpar[1].size());
+  Vector::iterator it=sValues.begin();
+  for (size_t j=0;j<gpar[1].size();++j) {
+    FiniteElement fe;
+    fe.v = gpar[1][j];
+    for (size_t i=0;i<gpar[0].size();++i) {
+      fe.u = gpar[0][i];
+      *it++ = field->valueFE(fe);
+    }
+  }
+  // Project the results onto the spline basis to find control point
+  // values based on the result values evaluated at the Greville points.
+  // Note that we here implicitly assume that the number of Greville points
+  // equals the number of control points such that we don't have to resize
+  // the result array. Think that is always the case, but beware if trying
+  // other projection schemes later.
+
+  RealArray weights;
+  if (surf->rational())
+    surf->getWeights(weights);
+
+  Go::SplineSurface* surf_new =
+        VariationDiminishingSplineApproximation(surf->basis(0),
+                                                surf->basis(1),
+                                                gpar[0], gpar[1],
+                                                sValues,
+                                                1,
+                                                surf->rational(),
+                                                weights);
+
+  vec.resize(surf_new->coefs_end()-surf_new->coefs_begin());
+  std::copy(surf_new->coefs_begin(),surf_new->coefs_end(),vec.begin());
+  delete surf_new;
+
+  return true;
+}
+
+
 // L2-Projection: Least-square approximation; global approximation
 Go::SplineSurface* ASMs2D::projectSolutionLeastSquare (const IntegrandBase& integrand) const
 {

src/ASM/ASMs3D.h

@@ -408,6 +408,11 @@ public:
   virtual bool evaluate(const ASMbase* basis, const Vector& locVec,
                         Vector& vec) const;
 
+  //! \copydoc ASMbase::evaluate(const Field*,Vector&)
+  //! \details Note a VDSA is used as the regular interpolation method in
+  //!          GoTools only works with uniform knots.
+  virtual bool evaluate(const Field* field, Vector& vec) const;
+
   //! \brief Evaluates the secondary solution field at all visualization points.
   //! \param[out] sField Solution field
   //! \param[in] integrand Object with problem-specific data and methods

src/ASM/ASMs3Drecovery.C

@@ -15,6 +15,8 @@
 #include "GoTools/trivariate/VolumeInterpolator.h"
 
 #include "ASMs3D.h"
+#include "FiniteElement.h"
+#include "Field.h"
 #include "CoordinateMapping.h"
 #include "GaussQuadrature.h"
 #include "SparseMatrix.h"
@@ -283,6 +285,60 @@ bool ASMs3D::globalL2projection (Matrix& sField,
 
 #include "ASMs3DInterpolate.C" // TODO: inline these methods instead...
 
+bool ASMs3D::evaluate (const Field* field, Vector& vec) const
+{
+  // Compute parameter values of the result sampling points (Greville points)
+  RealArray gpar[3];
+  for (int dir = 0; dir < 3; dir++)
+    if (!this->getGrevilleParameters(gpar[dir],dir))
+      return false;
+
+  // Evaluate the result field at all sampling points.
+  // Note: it is here assumed that *basis and *this have spline bases
+  // defined over the same parameter domain.
+  Vector sValues(gpar[0].size()*gpar[1].size());
+  Vector::iterator it=sValues.begin();
+  for (size_t l=0;l<gpar[2].size();++l) {
+    FiniteElement fe;
+    fe.w = gpar[2][l];
+    for (size_t j=0;j<gpar[1].size();++j) {
+      fe.v = gpar[1][j];
+      for (size_t i=0;i<gpar[0].size();++i) {
+        fe.u = gpar[0][i];
+        *it++ = field->valueFE(fe);
+      }
+    }
+  }
+
+  // Project the results onto the spline basis to find control point
+  // values based on the result values evaluated at the Greville points.
+  // Note that we here implicitly assume that the number of Greville points
+  // equals the number of control points such that we don't have to resize
+  // the result array. Think that is always the case, but beware if trying
+  // other projection schemes later.
+
+  RealArray weights;
+  if (svol->rational())
+    svol->getWeights(weights);
+
+  Go::SplineVolume* vol_new =
+        VariationDiminishingSplineApproximation(svol->basis(0),
+                                                svol->basis(1),
+                                                svol->basis(2),
+                                                gpar[0], gpar[1], gpar[2],
+                                                sValues,
+                                                1,
+                                                svol->rational(),
+                                                weights);
+
+  vec.resize(vol_new->coefs_end()-vol_new->coefs_begin());
+  std::copy(vol_new->coefs_begin(),vol_new->coefs_end(),vec.begin());
+  delete vol_new;
+
+  return true;
+}
+
+
 Go::SplineVolume* ASMs3D::projectSolutionLeastSquare (const IntegrandBase& integrand) const
 {
   if (!svol) return NULL;