Merge pull request #21 from pbauman/femfunction

Add System::project_solution/System::project_vector methods utilizing FEMFunctionBase functors

examples/adjoints/adjoints_ex3/coupled_system.h

@@ -114,7 +114,12 @@ class CoupledFEMFunctionsx : public FEMFunctionBase<Number>
   // Destructor
   virtual ~CoupledFEMFunctionsx () {}
 
- protected:
+  virtual AutoPtr<FEMFunctionBase<Number> > clone () const
+  {return AutoPtr<FEMFunctionBase<Number> >( new CoupledFEMFunctionsx(*this) ); }
+
+  virtual void operator() (const FEMContext&, const Point&,
+			   const Real, DenseVector<Number>&)
+  {libmesh_error();}
 
   virtual Number operator() (const FEMContext&, const Point& p,
 			     const Real time = 0.);
@@ -135,8 +140,13 @@ class CoupledFEMFunctionsy : public FEMFunctionBase<Number>
   // Destructor
   virtual ~CoupledFEMFunctionsy () {}
 
-
- protected:
+  virtual AutoPtr<FEMFunctionBase<Number> > clone () const
+  {return AutoPtr<FEMFunctionBase<Number> >( new CoupledFEMFunctionsy(*this) ); }
+
+  virtual void operator() (const FEMContext&, const Point&,
+			   const Real,
+			   DenseVector<Number>&)
+  {libmesh_error();}
 
   virtual Number operator() (const FEMContext&, const Point& p,
 			     const Real time = 0.);

include/numerics/const_fem_function.h

@@ -30,15 +30,24 @@ template <typename Output=Number>
 class ConstFEMFunction : public FEMFunctionBase<Output>
 {
 public:
-ConstFEMFunction (const Output c) {_c = c;}
+  ConstFEMFunction (const Output c) {_c = c;}
 
-~ConstFEMFunction() {}
+  ~ConstFEMFunction() {}
+
+  virtual AutoPtr<FEMFunctionBase<Output> > clone () const
+  {return AutoPtr<FEMFunctionBase<Output> >( new ConstFEMFunction(*this) ); }
 
   virtual Output operator() (const FEMContext&, const Point&,
 			     const Real /* time */ = 0.)
     { return _c; }
 
-  private:
+  virtual void operator() (const FEMContext&, const Point&,
+			   const Real,
+			   DenseVector<Output>& output)
+  {for(unsigned int i = 0; i < output.size(); i++ )
+      output(i) = _c;}
+
+private:
   Output _c;
 };
 

include/numerics/fem_function_base.h

@@ -47,16 +47,16 @@ class FEMFunctionBase
 protected:
 
   /**
-   * Constructor.  Optionally takes a master.
+   * Constructor.
    */
-FEMFunctionBase () {};
+  FEMFunctionBase () {};
 
 public:
 
   /**
    * Destructor.
    */
-virtual ~FEMFunctionBase () {};
+  virtual ~FEMFunctionBase () {};
 
 
   /**
@@ -65,8 +65,14 @@ virtual ~FEMFunctionBase () {};
    * Most problems will want to reimplement this for efficiency, in
    * order to call FE::get_*() as their particular function requires.
    */
-virtual void init_context (const FEMContext &) {}
+  virtual void init_context (const FEMContext &) {}
 
+  /**
+   * Returns a new copy of the function.  The new copy should be as
+   * ``deep'' as necessary to allow independent destruction and
+   * simultaneous evaluations of the copies in different threads.
+   */
+  virtual AutoPtr<FEMFunctionBase<Output> > clone () const = 0;
 
   // ------------------------------------------------------
   // misc
@@ -76,22 +82,74 @@ virtual void init_context (const FEMContext &) {}
    * Purely virtual, so you have to overload it.
    * Note that this cannot be a const method, check \p MeshFunction.
    */
-virtual Output operator() (const FEMContext&, const Point& p,
+  virtual Output operator() (const FEMContext&, const Point& p,
 			     const Real time = 0.) = 0;
 
 
+  /**
+   * Return function for vectors.
+   * Returns in \p output the values of the data at the
+   * coordinate \p p.
+   */
+  void operator() (const FEMContext&, const Point& p,
+		   DenseVector<Output>& output);
 
-protected:
+  /**
+   * Return function for vectors.
+   * Returns in \p output the values of the data at the
+   * coordinate \p p and for time \p time.
+   * Purely virtual, so you have to overload it.
+   * Note that this cannot be a const method, check \p MeshFunction.
+   */
+  virtual void operator() (const FEMContext&, const Point& p,
+			   const Real time,
+			   DenseVector<Output>& output) = 0;
+
+  /**
+   * @returns the vector component \p i at coordinate
+   * \p p and time \p time.
+   * Subclasses aren't required to overload this, since the default
+   * implementation is based on the full vector evaluation, which is
+   * often correct.
+   * Subclasses are recommended to overload this, since the default
+   * implementation is based on a vector evaluation, which is usually
+   * unnecessarily inefficient.
+   */
+  virtual Output component(const FEMContext&, unsigned int i,
+                           const Point& p,
+			   Real time=0.);
 
-/**
- * Variable index to decide which overloaded function should
- * be accessed
- */
 
-unsigned int var_index ;
+protected:
+
+  /**
+   * Variable index to decide which overloaded function should
+   * be accessed
+   */
+  unsigned int var_index ;
 
 };
 
+template <typename Output>
+inline
+Output FEMFunctionBase<Output>::component (const FEMContext& context, unsigned int i,
+					   const Point& p,
+					   Real time)
+{
+  DenseVector<Output> out(i+1);
+  (*this)(context, p, time, out);
+  return out(i);
+}
+
+template <typename Output>
+inline
+void FEMFunctionBase<Output>::operator() (const FEMContext& context, const Point& p,
+					  DenseVector<Output>& output)
+{
+  // Call the time-dependent function with t=0.
+  this->operator()(context, p, 0., output);
+}
+
 } // namespace libMesh
 
 #endif // LIBMESH_FEM_FUNCTION_BASE_H

include/systems/system.h

@@ -33,6 +33,7 @@
 #include "libmesh/qoi_set.h"
 #include "libmesh/reference_counted_object.h"
 #include "libmesh/variable.h"
+#include "libmesh/fem_function_base.h"
 
 // C++ includes
 #include <cstddef>
@@ -468,6 +469,18 @@ class System : public ReferenceCountedObject<System>
   void project_solution (FunctionBase<Number> *f,
                          FunctionBase<Gradient> *g = NULL) const;
 
+  /**
+   * Projects arbitrary functions onto the current solution.
+   * The function value \p f and its gradient \p g are
+   * user-provided cloneable functors.
+   * A gradient \p g is only required/used for projecting onto finite
+   * element spaces with continuous derivatives.
+   * If non-default \p Parameters are to be used, they can be provided
+   * in the \p parameters argument.
+   */
+  void project_solution (FEMFunctionBase<Number> *f,
+                         FEMFunctionBase<Gradient> *g = NULL) const;
+
   /**
    * Projects arbitrary functions onto the current solution.
    * The function value \p fptr and its gradient \p gptr are
@@ -499,6 +512,20 @@ class System : public ReferenceCountedObject<System>
                        FunctionBase<Number> *f,
                        FunctionBase<Gradient> *g = NULL) const;
 
+  /**
+   * Projects arbitrary functions onto a vector of degree of freedom
+   * values for the current system.
+   * The function value \p f and its gradient \p g are
+   * user-provided cloneable functors.
+   * A gradient \p g is only required/used for projecting onto finite
+   * element spaces with continuous derivatives.
+   * If non-default \p Parameters are to be used, they can be provided
+   * in the \p parameters argument.
+   */
+  void project_vector (NumericVector<Number>& new_vector,
+                       FEMFunctionBase<Number> *f,
+                       FEMFunctionBase<Gradient> *g = NULL) const;
+
   /**
    * Projects arbitrary functions onto a vector of degree of freedom
    * values for the current system.