NAME
FieldFE - finite element scalar field
INCLUDE
include "FieldFE.h"
SYNTAX
class FieldFE : public FieldWithPtValues
{
friend class FieldsFE;
protected:
Handle(GridFE) mesh;
Handle(Vec(NUMT)) nodal_values;
Handle(FiniteElement) globfe;
Handle(BasisFuncGrid) bfnodes;
Handle(ArrayGen(NUMT)) lattice_indexing; // same data as nodal_values
bool internal_basis_func_grid;
// precomputed data for improved efficiency of valueFEM/derivativeFEM:
int init4e; // the element for which FieldFE is precomputed
Vec(NUMT) loc_nodal_values; // element values of the field
bool different_bfgrid; // same bfgrid in FieldFE and FiniteElement?
bool different_geomgrid; // same grid in FieldFE and FiniteElement?
int nbf; // no of basis func./field d.o.f in an element
void precompute4interpolation (const FiniteElement& fe);
int ncalls_valuePt; // check that valuePt is not overused
bool reallocate
(
const GridFE& grid,
const Vec(NUMT)* coeff,
const BasisFuncGrid* sfg,
const char* fieldname
);
public:
FieldFE (); // needed by Field_prm::create function
FieldFE (const GridFE& grid, const char* fieldname);
FieldFE (const GridFE& grid, const Vec(NUMT)& coefficients,
const char* fieldname);
FieldFE (const BasisFuncGrid& bfnodes_grid, const char* fieldname);
FieldFE (const BasisFuncGrid& bfnodes_grid, const Vec(NUMT)& coefficients,
const char* fieldname);
~FieldFE ();
bool update (); // updating/redim due to a mesh change
bool compatible (const GridFE& g);
bool redim (const FieldFE& f, const char* fieldname);
bool redim (const GridFE& grid, const char* fieldname);
bool redim (const BasisFuncGrid& bfnodes_grid,
const char* fieldname);
bool redim (const GridFE& grid, const Vec(NUMT)& coefficients,
const char* fieldname);
bool redim (const BasisFuncGrid& bfnodes_grid,
const Vec(NUMT)& coefficients, const char* fieldname);
GridFE& grid () { return mesh.getRef(); }
const GridFE& grid () const { return mesh.getRef(); }
BasisFuncGrid& grid4basisFunc () { return bfnodes.getRef(); }
const BasisFuncGrid& grid4basisFunc () const { return bfnodes.getRef(); }
Vec(NUMT)& values () { return nodal_values.getRef();}
const Vec(NUMT)& values () const { return nodal_values.getRef();}
void localValues (VecSimple(NUMT)& local_values, int elm_no) const;
void fill (const Vector(NUMT)& new_values);
int getNoNodes () const;
bool isLattice (int nsd) const;
const ArrayGen(NUMT)& latticeIndex () const { return lattice_indexing(); }
ArrayGen(NUMT)& latticeIndex () { return lattice_indexing(); }
void exchangeValues (FieldFE& f); // this <-> f
virtual bool ok () const;
bool empty () const;
virtual void minmax (NUMT& min, NUMT& max, GridWithPts* grid =NULL) const;
NUMT valueElm (int elm_no, const Ptv(real)& local_pt, real t = DUMMY);
NUMT valueElm (int elm_no, const Ptv(real)& local_pt, real t = DUMMY) const
{ return CAST_CONST_AWAY(FieldFE)->valueElm(elm_no, local_pt,t); }
virtual NUMT valueFEM (const FiniteElement& fe, real t = DUMMY);
virtual NUMT valuePt (const Ptv(real)& x, real t = DUMMY);
NUMT valuePt (int& element, Ptv(real)& loc_pt, int& node,
const Ptv(real)& x, real t = DUMMY);
NUMT valuePt (const Ptv(real)& x, ElmSearchMethod k, real t);
NUMT valuePt ( int& element,
Ptv(real)& loc_pt,
int& node,
const Ptv(real)& x,
ElmSearchMethod k,
real t,
bool force_element_search = false,
int* n_Newton_iter = NULL,
real* Newton_error = NULL);
NUMT& valueNode (int node, real t = DUMMY); // inline
NUMT valueNode (int node, real t = DUMMY) const; // inline
virtual Ptv(NUMT) derivativePt (const Ptv(real)& x, real t = DUMMY);
Ptv(NUMT) derivativePt (const Ptv(real)& x, real t = DUMMY) const
{ return CAST_CONST_AWAY(FieldFE)->derivativePt(x,t); }
void derivativeNode (Ptv(NUMT)& gradient, int node, real t = DUMMY) const;
virtual void derivativeFEM
(Ptv(NUMT)& gradient, const FiniteElement& fe, real t = DUMMY);
void derivativeFEM
(Ptv(NUMT)& gradient, const FiniteElement& fe, real t = DUMMY) const
{ CAST_CONST_AWAY(FieldFE)->derivativeFEM (gradient, fe, t); }
void derivativeElm (Ptv(NUMT)& gradient, int elm_no,
const Ptv(real)& local_pt, real t = DUMMY);
void derivativeElm (Ptv(NUMT)& gradient, int elm_no,
const Ptv(real)& local_pt, real t = DUMMY) const
{ CAST_CONST_AWAY(FieldFE)->derivativeElm (gradient, elm_no, local_pt, t); }
// interpolate: op= for fields of different size (over different grids)
void interpolate (const FieldFE& fefield);
void interpolate (const FieldLattice& fdfield);
// op= implies a common underlying grid for the two fields
// (no interpolation is needed, i.e., op= is much faster than interpolate)
void operator = (const FieldFE& fefield);
void operator = (const FieldLattice& fdfield);
void operator = (const FieldFunc& func);
void fill (const FieldFunc& func, real time);
virtual void fill (NUMT value);
virtual void add (NUMT value);
virtual void mult (NUMT value);
virtual void apply (Func(NUMT) f);
virtual void add (Field& field, int power, NUMT front_factor);
void add (const FieldFE& f); // more efficient this += f
virtual Field& scale ();
virtual Field& unscale ();
virtual void unloadData (Os os) const;
virtual void loadData (Is is);
virtual void attach (Grid& grid);
virtual Grid* getGridBase();
// virtual functions from FieldWithPtValues:
virtual int getNoPoints () const;
virtual int getNoValues () const;
virtual NUMT& valuePoint (int point_no)
{ return nodal_values()(point_no); } // never inline because of virtual
virtual NUMT valuePoint (int point_no) const
{ return nodal_values()(point_no); }
virtual Ptv(real) getPt (int point_no) const
{ return bfnodes->getCoor(point_no); }
virtual GridWithPts& getGridWithPts ()
{ return mesh.getRef(); }
virtual const GridWithPts& getGridWithPts () const
{ return mesh.getRef(); }
virtual Vec(NUMT)& valuesVec ()
{ return nodal_values.getRef(); }
virtual const Vec(NUMT)& valuesVec () const
{ return nodal_values.getRef(); }
CLASS_INFO
VIRTUAL_CAST(FieldFE);
// temporary:
COPY_CONSTRUCTOR(FieldFE);
virtual void print (Os os) const;
void scan (Is is);
};
KEYWORDS
finite element, scalar field, FieldFE
DESCRIPTION
The class implements a scalar field defined in terms of a finite
element grid and assocated finite element interpolation func
tions. That is, the field consists of a GridFE object and the
coefficients in the finite element expansion. The class is funda
mental when coding simulators based on finite elements in Diff
pack.
CONSTRUCTORS AND INITIALIZATION
One can call the constructor without arguments and afterwards
initialize the object by calling a relevant redim function. As
usual in Diffpack, "initialization" here refers to allocating
dynamic memory and not to filling the nodal values.
If the grid consists of isoparametric elements only, one can use
the constructors that take a GridFE object. If the vector of
nodal values is available (for example, if several fields shall
share the same vector) it can be given as an optional argument to
the constructor. For element meshes containing non-isoparametric
elements, one must use the constructor (or redim function) that
takes a BasisFuncGrid as argument. Each constructor has a corre
sponding redim function.
MEMBER FUNCTIONS
update - redimensions the nodal values array (and the internal
BasisFuncGrid object if that object was made by this field - oth
erwise the simulator or another field that made the BasisFuncGrid
object will update it) if the grid has changed, for example, if
the number of nodes is changed. It is important to call this
function when one uses adaptive grids.
compatible - checks if this field is compatible with the grid
given as argument, that is, if this field's grid is exactly the
same as the argument grid and if this field's nodal array has the
same size as the number of nodes in the argument grid.
grid - returns access to the underlying finite element grid.
grid4basisFunc - returns access to the information on basis func
tion nodes. In
case of isoparametric elements one can use "grid" instead,
but
a call to grid4basisFunc is more general since it can handle the
case where nodes defining the geometry and the basis functions
differ. The BasisFuncGrid object is convenient when working with
mixed finite elements.
values - returns access to the vector Vec(NUMT) (that is,
Vec(real) or Vec(Complex)) of the degrees of freedom in the
field. For isoparametric elements this vector contains the nodal
values of the field. For example, if one wants to print the nodal
values, one can then say
field.values().print(cout,"nodal values of field")
(The FieldFE::print function prints both the grid and the nodal
values, so the construction above is convenient to print the val
ues only.) For more general finite element fields over non-
isoparametric elements where there are more than one degree of
freedom at the nodes, one must use the function valueNode to
extract the values at a node, since some of the degrees of free
dom may not be interpreted as nodal values.
localValues - extracts the entries in nodal_values (alterna
tively, the entries in the vector returned by values()) that are
associated with an element elm_no. local_values(i) is then, for
isoparametric elements, the field value at local node i in ele
ment elm_no. Another way of extracting this value is f.val
ues()(f.grid4basisFunc().loc2glob(elm_no,i)) for a FieldFE f.
(However, it is usually more efficient first to extract
local_values and then index this short vector.)
fill - fills the nodal values vector with values given by the
vector argument.
getNoNodes - returns the number of basis function nodes. Same as
getNoPoints.
isLattice - returns true if the underlying finite element mesh is
a regular lattice with uniform spacing. Then it is possible to
use or indexing of nodes by calling latticeIndex. The function
takes and int nsd argument. The reason for this is that if the
grid is a lattice, the user must know the number of space dimen
sions in order to index the nodal values vector (which in that
case is an ArrayGen(real) object) correctly. The argument checks
that the user knows enough to use the lattice grid facility.
latticeIndex - returns access to an ArrayGen(real) array contain
ing the nodal values. This enables "finite difference" indexing
of the nodes in a finite element grid. There is no extra storage
associated with this functionality (values and latticeIndex
return access to the same physical array). Always call isLattice
prior to latticeIndex to ensure that the finite element grid
really is a lattice. If isLattice returns false, latticeIndex
will return reference to a NULL pointer. For efficiency reasons,
there is no warning message in this case.
minmax - finds the minimum and maximum nodal values in the field.
valueElm - interpolates the field value at a local point in an
element.
valueFEM - interpolates the field value at a local point in an
element. (valueFEM is the recommended function, valueElm is only
included for backward compatibility).
valuePt - computes the value of the finite element field at a
spatial point. The time variable has no effect in the current
implementation. There are several overloaded versions of val
uePt, each offering different search algorithms. In fact, GridFE
functions search for the point and in FieldFE we only carry out
the interpolation. Using +verbose n on the command line, where n
is an integer greater than zero, most search algorithms will give
messages to the standard output if the search was unsuccessful.
This is an important feature if you experience strange numerical
results in a code that calls valuePt. An overloaded version of
valuePt enables the user to choose between several search algo
rithms for finding the element and the local point. See the
function GridFE::findElmAndLocPt for more information about the
various choices.
valueNode - gives access to the finite element field at a node.
Notice that this function can also be used on the left hand side
of an equality operator (for assigning nodal values in a field).
For isoparametric elements one can also call values and then
invoke the operator() function of the returned Vec(NUMT) object.
Note that valueNode is not inline because it needs to check if
there are more than one field degree of freedom at each node.
Therefore, valueNode(i) is very much less efficient than val
ues()(i), but the latter may be wrong if one works with certain
types of non-isoparametric elements (more than one unknown per
node).
derivativePt - as valuePt but the function computes the deriva
tive.
derivativeNode - the counterpart to valueNode, but when the ele
ments have discontinuous derivatives at the element boundaries,
this function is not very useful. It is included just to give the
user an error message about what should be done to compute the
derivative at a node. The discontinuous derivative field should
first be computed and then smoothed and represented as a new
FieldFE object which can be evaluated at a node by calling val
ueNode.
derivativeFEM - as valueFEM, but the gradient is computed.
derivativeElm - as derivativeFEM but other arguments are offered
(the element number and the local point in that element).
interpolate - given another field on an arbitrary grid, this
function computes the nodal values of the "this" grid by standard
interpolation (the function valuePt is used). Note that "this"
refers to the f field in the call f.interpolate(g). The interpo
late function requires the "this" field to be initialized with a
grid. If you have an empty field (without an associated grid)
and want to initialize it with another field, use the operator=
functions. operator= is intended for setting two fields equal to
each other, where also the grids are identical. interpolate can
set two fieldss equal to each other where the grids differ. Of
course, operator= is very much more efficient than interpolate
(since it avoids the use of valuePt and instead can only copy
array entries).
operator=(FieldLattice) - enables to convert a finite difference
field FieldLattice to a finite element field. Each cell in the
difference grid becomes a ElmTensorProd1 element.
operator=(FieldFunc) - fills the nodal values by evaluating an
explicit function (either a standard function or a C++ functor is
offered by class FieldFunc).
fill(FieldFunc, real ) - as operator= above, but a time parameter
can be given. This is useful when the explicit function changes
with time. The operator= function goes through all nodes and
calls the explicit function with a DUMMY argument for the time
variable (since operator= can only take one argument). fill is
intended for the situations where one needs two arguements, which
is the general case. It is a good habit to use the fill function
instead of the fancy assignment syntax when filling finite ele
ment fields with explicit function values. This makes the exten
sion to time dependent functions easier.
ok - returns a true value if the object is in an ok state, that
is, if all the internal data structures are allocated. Otherwise,
it gives some error messages.
empty - returns true if the object has not allocated the internal
data structures. In contrast to ok, no error messages are issued.
Hence the function is convenient in tests where one wants to call
a redim function if the field is empty.
The rest of the functions are virtual and explained in the base
classes.
SEE ALSO
class Field, class GridFE
DEVELOPED BY
SINTEF Applied Mathematics, Oslo, Norway, and University of Oslo,
Dept. of Mathematics, Norway
AUTHOR
Hans Petter Langtangen, SINTEF/UiO