// ===========================================================================
//
// Filename: FEM4EllipticPDE.h
//
// Created: 08/16/12 18:19:41
// Compiler: Tested with g++, icpc, and MSVC 2010
//
// Author: Trevor Irons (ti)
//
// Organisation: Colorado School of Mines (CSM)
// United States Geological Survey (USGS)
//
// Email: tirons@mines.edu, tirons@usgs.gov
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// ===========================================================================
/**
@file
@author Trevor Irons
@date 08/16/12
@version $Rev$
**/
#ifndef FEM4ELLIPTICPDE_INC
#define FEM4ELLIPTICPDE_INC
#include "LemmaObject.h"
#include "vtkImplicitFunction.h"
#include "vtkDataSet.h"
#include "vtkXMLDataSetWriter.h"
#include "vtkSmartPointer.h"
#include "vtkIdList.h"
#include "vtkCell.h"
#include "vtkDoubleArray.h"
#include "vtkPointData.h"
#include "vtkDataSetSurfaceFilter.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include
#include "DCSurvey.h"
namespace Lemma {
/** @defgroup FEM4EllipticPDE
@brief Finite element solver for elliptic PDE's
@details Uses Galerkin method.
*/
/** \addtogroup FEM4EllipticPDE
@{
*/
// ===================================================================
// Class: FEM4EllipticPDE
/**
@class FEM4EllipticPDE
@brief Galerkin FEM solver for Elliptic PDE problems
@details Solves problems of the type
\f[ - \nabla \cdot ( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) ) = g(\mathbf{r})
\f]
*/
class FEM4EllipticPDE : public LemmaObject {
friend std::ostream &operator<<(std::ostream &stream,
const FEM4EllipticPDE &ob);
public:
// ==================== LIFECYCLE =======================
/** Returns a pointer to a new object of type FEM4EllipticPDE.
* It allocates all necessary memory.
*/
static FEM4EllipticPDE* New();
/** Deletes this object. Delete also disconnects any
* attachments to this object.
*/
void Delete();
// ==================== OPERATORS =======================
// ==================== OPERATIONS =======================
/** Sets the vtkImplictFunction that is used as sigma in
\f$ - \nabla \cdot ( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) ) = g(\mathbf{r})
\f$. If this is not set, sigma evaluates as 1.
*/
void SetSigmaFunction(vtkImplicitFunction* sigma);
/** Sets the step to be used in applying boundary conditions.
*/
void SetBoundaryStep(const Real& h);
/** Sets the vtkImplictFunction that is used
*/
void SetGFunction(vtkImplicitFunction* G);
/** Uses a function pointer instead of a vtkImplicitFunction. For functions that can
* be written in closed form, this can be much faster and more accurate then interpolating
* off a grid.
*/
void SetGFunction( Real (*pFcn3)(const Real&, const Real&, const Real&) );
/** Sets the vtkDataSet that defines the calculation grid.
*/
void SetGrid(vtkDataSet* Grid);
/** Sets up a DC problem with a survey
* @param[in] ij is the injection index
*/
void SetupDC(DCSurvey* Survey, const int& ij);
/** Performs solution */
void Solve(const std::string& fname);
/** Performs solution */
void SolveOLD(const std::string& fname);
// ==================== ACCESS =======================
// ==================== INQUIRY =======================
protected:
// ==================== LIFECYCLE =======================
/// Default protected constructor.
FEM4EllipticPDE (const std::string& cname);
/** Default protected constructor. */
~FEM4EllipticPDE ();
/**
* @copybrief LemmaObject::Release()
* @copydetails LemmaObject::Release()
*/
void Release();
// ==================== OPERATIONS =======================
/** Used internally to construct stiffness matrix, A
*/
void ConstructAMatrix();
/* Used internally to construct the load vector, g
*/
//void ConstructLoadVector();
// This function is copyright (C) 2012 Joseph Capriotti
/**
Returns a vtkIdList of points in member vtkGrid connected to point ido.
@param[in] ido is the point of interest
@param[in] grid is a pointer to a vtkDataSet
*/
vtkSmartPointer GetConnectedPoints(const int& id0);
/** Uses Simpson's rule to numerically integrate
* (Shamelessly adapted from http://en.wikipedia.org/wiki/Simpson's_rule
*/
Real CompositeSimpsons(vtkImplicitFunction* f, Real r0[3], Real r1[3], int numInt);
/** Uses Simpson's rule to numerically integrate
* (Shamelessly adapted from http://en.wikipedia.org/wiki/Simpson's_rule
* this method is for a static f
*/
Real CompositeSimpsons(const Real& f, Real r0[3], Real r1[3], int numInt);
/**
* Uses Simpson's rule to numerically integrate two functions in 1 dimension over the points
* r0 and r1
*/
Real CompositeSimpsons2(vtkImplicitFunction* f, Real r0[3], Real r1[3], int numInt);
/**
* Uses Simpson's rule to numerically integrate two functions in 1 dimension over the points
* r0 and r1, uses the Hat function and the function pointer fPtr
*/
Real CompositeSimpsons2( Real (*fPtr)(const Real&, const Real&, const Real&), Real r0[3], Real r1[3], int numInt);
/**
* Uses Simpson's rule to numerically integrate two functions in 1 dimension over the points
* r0 and r1, uses the Hat function. Assumes a constand function value f
*/
Real CompositeSimpsons2( const Real& f, Real r0[3], Real r1[3], int numInt);
/**
* Uses Simpson's rule to numerically integrate two functions in 1 dimension over the points
* r0 and r1, uses the Hat function. Assumes a linear function from f0 to f1.
*/
Real CompositeSimpsons2( const Real& f0, const Real& f1, Real r0[3], Real r1[3], int numInt);
/**
* Uses Simpson's rule to numerically integrate three functions in 1 dimension over the points
* r0 and r1, uses two Hat functions. Assumes a linear function from f0 to f1.
*/
Real CompositeSimpsons3( const Real& f0, const Real& f1, Real r0[3], Real r1[3], int numInt);
/** Hat function for use in Galerkin FEM method, actually a 1/2 Hat.
*/
//Real hat(Real r[3], Real ri[3], Real rm1[3] );
Real Hat(const Vector3r &r, const Vector3r& ri, const Vector3r& rm1 );
/** Computes distance between r0 and r1
*/
Real dist(Real r0[3], Real r1[3]);
/** Computes distance between r0 and r1
*/
Real dist(const Vector3r& r0, const Vector3r& r1);
// ==================== DATA MEMBERS =========================
/** The effective step at the boundary condition. Defaults to 1
*/
Real BndryH;
/** The sigma value at the boundary condition. Defaults to 1
*/
Real BndrySigma;
/** Implicit function defining the \f$\sigma\f$ term in
\f$ - \nabla \cdot ( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) ) = g(\mathbf{r})
\f$.
If set to 1, the Poisson equation is solved. Any type of vtkImplicitFunction may
be used, including vtkImplicitDataSet.
*/
vtkSmartPointer vtkSigma;
/** Implicit function defining the \f$g\f$ term in
\f$ - \nabla \cdot ( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) ) = g(\mathbf{r})
\f$.
Typically this is used to define the source term.
*/
vtkSmartPointer vtkG;
/** Any type of vtkDataSet may be used as a calculation Grid */
vtkSmartPointer vtkGrid;
/** Function pointer to function describing g */
Real (*gFcn3)(const Real&, const Real&, const Real&);
Eigen::SparseMatrix A;
VectorXr g;
private:
}; // ----- end of class FEM4EllipticPDE -----
/*! @} */
} // ----- end of Lemma name -----
#endif // ----- #ifndef FEM4ELLIPTICPDE_INC -----