NAME
GridFE - finite element grid
INCLUDE
include "GridFE.h"
SYNTAX
enum ElmSearchMethod { // method for GridFE::findElmAndLocPt functionality
MAX_GRID_SPACING, // see GridFE.cpp for explanation of the enum values
LATTICE_GRID,
EBE_SEARCH
};
class GridFE : public GridWithPts
{
friend class FieldFE; // for faster interpolation
friend class BasisFuncGrid; // for faster interpolation
friend class FieldPiWisConst;
// some subclasses need to be friends:
friend class GridFEInfo;
friend class GridFEAdB;
friend class GridDynFE;
// some preprocessors need to be friends:
friend class PreproBox;
friend class PreproGeomPack;
friend class ManipGridDynFE;
protected:
MatSimple(int) nodel; // nodes in the elements
MatSimple(real) coor; // coordinates of the nodes
Handle(Indicators) bind; // boundary indicators and names
int nne; // number of nodes in an element
VecSimple(int) nne_all; // if nne=0, then use nne_all(1:nel) array
int maxnne;
String elm_tp; // type of element (geometry) (if nne>0)
VecSimple(String) elm_tp_all; // if nne=0, use elm_tp_all(nel) for elm_tp
bool onemat; // indicator: ONE or many materials
VecSimple(int) material_type; // material type array, used if onemat==0
int nsd; // number of element space dimensions
int nsd_phys; // grid embedded in nsd_phys space dim.
int nno; // number of nodes in the grid
int nel; // number of elements in the grid
public:
GridFEInfo additional_info; // additional grid data structures
BasisFuncGrid* associated_bfg; // may point to a b.f.g. binded to this grid
Handle(Indicators) elmind; // for marking elements for subdomains
protected:
bool uniform_mesh; // true: all elements are equal in shape
bool same_elm_tp; // true: all elements are of the same type
Handle(GridLattice) lattice_data;
HandleElmDefs elmdef;
public:
virtual int getNoSpaceDim () const;
int getNoPhysicalSpaceDim () const;
int getNoNodes () const { return nno; }
int getNoElms () const { return nel; }
int getNoBoInds () const { return bind->getNoInds(); }
bool allElmsOfSameType () const { return getbool(nne > 0); }
bool hasUniformMesh () const { return uniform_mesh; }
void setUniformMesh ();
void setNonUniformMesh ();
int getNoMaterials () const; // requires significant computations!
GridFE();
GridFE (const GridFE& grid);
bool redim
(
int nsd, // no of space dimensions
int nno, // no of nodes
int nel, // no of elements
int maxnne, // max no of nodes in an element
int nbind, // no of boundary indicators
int nne, // =0: variable no of nodes in elements, or =maxnne
bool onemat = true // one material (true) or several (false)
);
~GridFE();
bool ok () const;
void operator = (const GridFE& grid);
virtual bool notEqual (const GridFE& grid, int check_level = 2,
bool silent = false);
int getNoNodesInElm (int e) const {return nne==0 ? nne_all(e):nne; }
int getMaxNoNodesInElm () const { return maxnne; }
int getMaterialType (int e) const;
bool oneMaterialType() { return onemat; }
String getElmType (int e) const;
void getElmCoor (Mat(real)& elmcoor, int elm_no) const;
real getCoor (int node, int dir) const { return coor(node,dir); }
Ptv(real) getCoor (int node) const;
void getCoor (Ptv(real)& point, int node) const;
virtual void getMinMaxCoord (Ptv(real)& mincoord, Ptv(real)& maxcoord) const;
Ptv(real) getLocalCentroid (int elm_no);
Ptv(real) getGlobalCentroid (int elm_no);
Ptv(real) getGlobalMeanCentroid (int elm_no);
// local <-> global node numbering mapping:
int loc2glob (int e, int local_node) const { return nodel(e,local_node);}
#define loc2glob_eff loc2glob /* for backward compatibility */
// boundary indicator information:
String getBoIndName (int indno) const { return ( bind->getName(indno) ); }
// are all nodes in side side in element elm_no marked with
// boundary indicator no. indno ?
bool boSide (int elm_no, int side, int indno, const ElmDef& elm) const;
bool boNode (int node, int indno) const { return bind->on(node,indno); }
bool boNode (int node) const; // is node marked with any indicator?
// is any node in element elm_no marked with indicator indno?
bool boNodesInElm (int elm_no, int indno) const;
bool boNodesInElm (int elm_no) const; // marked with any indicator?
bool boInd (int ind) const; // any node with indicator no. ind?
#define BoSide boSide /* for backward compatibility */
#define BoNode boNode /* for backward compatibility */
#define BoNodesInElm boNodesInElm /* for backward compatibility */
#define BoInd boInd /* for backward compatibility */
const Indicators& boundaryData () const { return bind.getRef(); }
Indicators& boundaryData () { return bind.getRef(); }
void attach (const Indicators& boinds);
// modify boundary indicators:
void redefineBoInds (Is is) { bind->redefine (is); }
void addBoIndNodes (Is is, bool add = ON /*OFF: delete indicators*/);
void addMaterial (Is is);
int addMaterial(const int material_number,const MatSimple(real)& polygon);
int insidePolygon(const MatSimple(real)& polygon, const Ptv(real)& globpt );
// lattice functionality (uniform partition, identical box elements):
void scanLattice (Is is);
void operator= (const GridLattice& lattice);
virtual bool isLattice () const;
virtual bool isReallyLattice () const; // more thorough check than isLattice
const GridLattice& getLattice () const { return lattice_data(); }
void refineIfBox (const Ptv(int)& refinement_factor);
// users should only apply the following move functions to change
// the nodal coordinates (and not putCoor):
void move (const FieldsFE& displacement, real scale_factor = 1.0);
void move (FieldsFunc& func, real time = DUMMY);
void move (const Mat(real)& new_coor);
void move (const Vec(real)& new_coor, int dir);
void move (const GridFE& grid);
// iterate through all the grid points:
virtual void startIterator ();
virtual bool nextPt (Ptv(real)& x);
bool nextPt (Ptv(real)& x, int& node);
virtual int getNoPoints () const { return nno; }
// search for points and nodes in elements:
int isNode (const Ptv(real)& point);
virtual int nearestPoint
(const Ptv(real)& point,
real& distance,
bool& exact);
bool isNodeInElm (int node, int e);
void boundingBox (int e, Ptv(real)& min, Ptv(real)& max);
void findElmAndLocPt // offers different algorithms (k arg.)
(
const Ptv(real)& glob_pt,
int& element,
Ptv(real)& loc_pt,
int& node,
ElmSearchMethod k,
bool force_element_search = false,
int* n_Newton_iter = NULL,
real* Newton_error = NULL);
void findElmAndLocPt // original
(
const Ptv(real)& glob_pt,
int& element,
Ptv(real)& loc_pt,
int& node,
bool safe_procedure = false, // efficient, backw comp.
bool force_element_search = false,
int* n_Newton_iter = NULL,
real* Newton_error = NULL);
// input/output:
friend Os& operator << (Os& os, const GridFE& mesh);
friend Is& operator >> (Is& is, GridFE& mesh);
virtual void print (Os os) const;
virtual void scan (Is is);
virtual void scale (); // scale the nodal coordinates
virtual void unscale (); // unscale the nodal coordinates
void attachScale (const SpaceTimeScale& scale);
// useful when comparing grid points with a tolerance:
real calcTolerance (real div_factor = 100.0) const;
// -------------------------------------------------------------
// functions for setting grid data, mostly used by preprocessors
void setMaterialType (int e, int material_number);
void setElmType (int e, const String& type); // for heterogeneous grids
void setElmType (const String& from, const String& to);
void setElmType (const String& name); // if allElmsOfSameType
void setNoNodesInElm (int e, int nnodes);
void fillCoor (real value) { coor.fill(value); }
void setNoPhysicalSpaceDim (int nsd_phys);
// putCoor and putLoc2Glob should not be call by preprocessors, not users!
void putCoor (real value, int node, int dir) { coor(node,dir)=value; }
void putCoor (const Ptv(real)& point, int node);
void putLoc2Glob (int global_node, int e, int local_node)
{ nodel(e,local_node) = global_node; }
void setBoInd (int node, int indno) { bind->set(node,indno); }
void clearBoInd (int node, int indno) { bind->clear(node,indno); }
void clearBoInds () { bind->fill(OFF); }
void putBoIndName (const String& name, int indno)
{ bind->putName(name,indno); }
int searchPoint (const Ptv(real)& point);
void updateElmConnect( int nodeno, int newnodeno);
virtual void newNodalNumbering (VecSimple(int)& old2new);
virtual void isToBeChanged (); // update associated_bfg if necessary
// for adaptive grids:
void nodeCouplings (int node, VecSimple(int)& couplings);
virtual void constraintCouplings(VecSimplest(SetDistinct(int))& /*n2n*/,
int /*array_no*/ = 0) {}
virtual int getNoConstraints () const { return 0; }
virtual int calcBandwidth () const;
virtual int getGridLevelNo () { return 1; }
virtual GridFE& getStartGrid() { return *this; }
virtual void modify4constraints
(DegFreeFE& /*dof*/, Matrix(NUMT)& /*A*/, Vector(NUMT)& /*b*/) {}
virtual void modifyField4constraints (FieldFE&) {}
protected:
void printAscii (Os os) const;
void scanAscii (Is is);
void printBinary (Os os) const;
void scanBinary (Is is);
int current_node; // used in the grid iterator
Mat(real) elmcoor_; // just to avoid alloc. in some func.
public:
CLASS_INFO
VIRTUAL_CAST(GridFE)
};
KEYWORDS
finite element, grid, mesh, GridFE
DESCRIPTION
The class represents a finite element grid with information on
nodal coordinates, the element topology and the nodes that are
marked with boundary indicators.
Class GridFE objects are usually declared as empty objects and
filled by preprocessor classes. Finite element computations
require normally access to the read-only functions in class
GridFE. For preprocessor programmers and for educational pur
poses direct use of the assignment functions in the GridFE class
is necessary.
CONSTRUCTORS AND INITIALIZATION
There is an empty constructor (that does not allocate memory for
the internal data structure) and a constructor that takes some
parameters like the number of nodes and elements and allocates
the memory for the internal data structure. When using the empty
constructor, allocation of memory must be performed later by the
redim function. Notice that this function only allocates memory.
The coordinates, the element connectivity, the element types and
so must be initialized afterwards. To fill the data structures,
one can use the many member functions, like e.g. putCoor, or one
can use scan to read a file in a certain format. Examples on this
format are given in the examples below. One can also use the
scanLattice function that reads a string having GridLat
tic::scan(Is) syntax and initializes a lattice grid. There is
also an operator= (const GridLattice function available that can
be used to initialize the finite element grid if one has a uni
form finite difference at hand (remember that the GridLattice
object sent to operator= must be handled by a handle such that it
is not deleted - this grid object binds a handle to the GridLat
tice object in operator=).
BASIC FUNCTIONALITY
A finite element grid consists of nodes and elements. The number
of space dimensions, elements and nodes are given by the func
tions getNoSpaceDim, getNoNodes and getNoElms. The number of
nodes in an element is returned by getNoNodesInElm. The coordi
nates of the nodes can be read and set using the member functions
getCoor and putCoor. Notice that there are several overloaded
versions of these functions. The applications of the functions
will determine which version that may be most convenient to use.
If an existing grid is to be deformed, use the move functions
rather than the putCoor functions. The move functions are safer.
For example, if a node is moved by putCoor and this makes an ini
tially uniform mesh non-uniform, one must also call the set
NonUniformMesh function. If not, class FiniteElement will find
out that the grid is uniform and apply particularly efficient
techniques when computing finite element data. These data will
be completely wrong if the mesh is non-uniform. If NUMT is Com
plex, the move function corresponding to a displacement field
will only use the real part of that field for moving the nodal
points.
The nodes that belong to the various elements are found from the
member function loc2glob, or more precisely, given a local node
number in an element, loc2glob returns the corresponding global
number of that node. To set a new value of the global node number
one may use the function putLoc2Glob. The function getCentroid
returns the local element coordinates of the centroid of an ele
ment. Convenient in combination with class FieldPiWisConst.
It may be convenient to mark some nodes in a grid where boundary
conditions are to be applied, or the marking may have applica
tions when plotting the grid. The GridFE class therefore supports
marking of nodes: a node can be marked with a so called boundary
indicator. Using class Indicators, its use is transparent to
users of class GridFE, one can define a set of boundary indica
tors. Each indicator is given a name by using the putBoIndName
function. Then the member function setBoInd can be used to set a
boundary indicator in ON mode at a particular node. To turn on
the OFF mode, call clearBoInd. To test whether a node is marked
with a particular boundary indicator, one can call boNode which
returns true if the node is marked, that is, if the particular
boundary indicator is in ON mode. The function boNodes returns
true if any node in an element is marked with a particular bound
ary indicator. The BoInd function can be used to check if a par
ticular indicator is on in any node. For example, if an applica
tion does not allow the use of a particular boundary condition,
it is easy to check if the user has set an indicator for this
condition by using BoInd.
Normally, the GridFE object allocates an internal Indicators
object for keeping track of boundary indicators. The programmer
can at any time replace this object by another one using the
attach function. For example, if several simulators share the
same grid, but need different boundary indicators, the programmer
can have a vector of Indicators objects and attach these to the
grid object as necessary. See the FAQ for an example.
The domain can be partitioned into materials, which are identi
fied as distinct (non-overlapping) subdomains. These might, for
example, correspond to different physical media (different mate
rials in heat transfer and elasticity). Each element can hence
be given a material number. To get and set the material number
of an element the member functions getMaterialType and setMateri
alType are used. If there is only one material it is not neces
sary to set the subdomain type and the internal data structure of
class GridFE can be made simpler. The function oneMaterialType is
used to see if only one material type is used. The redim function
requires an argument to indicate whether there are multiple mate
rials. Even if the grid is initialized as one material, one can
easily add new material subdomains using the addMaterial func
tions.
It is emphasized that "material" and "subdomain" are two differ
ent concepts in class GridFE. While "material" refer to non-over
lapping parts of the domain, with different physical materials,
the "subdomain" concept applies to possibly overlapping subdo
mains, which are useful for domain decomposition methods and par
allel computing. Subdomains are represented by the elmind indi
cators, whereas materials are represented by the material_type
array. The rule of thumb is: One element can belong to several
subdomains, but only one material.
In GridFE one associates an element geometry type with each ele
ment. The element geometry type is determined from the functions
that are used in the mapping from the local element coordinate
system onto the physical coordinate system. The functions setElm
Type and getElmType are used to set and read the geometry type.
The variable in this context is a string that coincides with the
name of the element. This name is identical to the classname of
a subclass in the ElmDef hierarchy. In Diffpack the definitions
of various finite elements are given as objects in the ElmDef
hiearchy. Examples of classes derived from ElmDef are
ElmTensorProd1: Multi-dimensional tensor product elements, where
the underlying one-dimensional basis functions are linear (the
last number, 1, indicates first order). The classname must be
supplied by the number of space dimensions in order to specify
the element completely.
ElmTensorProd2: As ElmTensorProd1, but the underlying one-
dimensional basis functions are quadratic.
ElmTensorProd: As ElmTensorProd1 and ElmTensorProd2, but the
order of the underlying one-dimensional basis function is a
required parameter for the object (in addition to the number of
space dimension).
ElmB2n1D: Element of box shape, with 2 nodes, in 1D. This is the
same element as ElmTensorProd1 in 1D or ElmTensorProd in 1D with
linear underlying one-dimensional basis functions. However, the
implementation is much more efficient; in 1D cases ElmTensorProd
sets up a lot of loops of unit length, while ElmB2n1D computes
the two basis functions straightforwardly. In CPU intensive
applications one should therefore avoid the general, but some
times very convenient, ElmTensorProd family of elements and
instead use the specialized implementations such as ElmB2n1D.
ElmB4n2D: Another specialization of ElmTensorProd or ElmTensor
Prod1: Element of box shape, with 4 nodes, in 2D - that is, the
standard bilinear quadrilateral element. The implementation of
ElmB4n2D is much more efficient than the general class ElmTensor
Prod (or ElmTensorProd1), see the comments under the ElmB2n1D
element.
ElmB9n2D: Another specialization of ElmTensorProd or ElmTensor
Prod2: Element of box shape, with 9 nodes, in 2D.
ElmB8n2D: Eight node, quadrilateral element (as ElmB9n2D, but
without the internal node).
ElmB8n3D: A specialization of ElmTensorProd and ElmTensorProd1
in 3D: The standard 8 node trilinear brick element. Again this
implementation is more efficient than the ElmTensorProd code, see
comments under ElmB2n1D.
ElmT3n2D: Element of triangular shape, with 3 nodes, in 2D.
ElmT4n3D: The standard tetrahedra (linear) element with 4 nodes
in 3D.
For example, if the element is of bilinear type, getElmType will
return "ElmB4n2D" (or the equivalent "ElmTensorProd", depending
on what was set using the setElmType function). Consult the
ElmDef hierarchy to see what types of elements that are sup
ported. In many cases all elements are of the same type, and
setElmType needs to be called only once.
If all elements are isoparametric, the geometry element type is
sufficient in all instances for determining the type of finite
element. If there are non-isoparametric elements in the mesh,
class BasisFuncGrid defines the nodes for the trial functions and
overrides many of the functions in class GridFE. In particular,
the element type is redefined in BasisFuncGrid if there are non-
isoparametric elements.
MEMBER FUNCTIONS
Many member functions are explained above. A short description is
given below.
redim - sets new dimensions of the internal data structure in the
object. It is necessary to call redim if the number of nodes or
elements has changed. All old values in the grid object will be
forgotten after a call to redim so a complete new initialization
is required.
getNoNodesInElm - get the number of nodes in an element.
setNoNodesInElm - set the number of nodes in an element.
getMaterialType - get the material number of an element.
setMaterialType - set the material number of an element.
getElmType - get the element geometry type (a string is returned
containing the name of the element, where the name is the name of
a subclass in the ElmDef hierarchy.
setElmType - set the element geometry type. Can be called after
the grid is generated by a preprocessor to change the element
type to, e.g., a new element made particularly for an applica
tion.
fillCoor - initializes the coordinates of all nodes (usually this
function is used to set all coordinate values to zero prior to
the computation of a grid).
getElmCoor - get the coordinates of all nodes in an element.
getCoor - get the coordinates of a node. There are several over
loaded versions. Maximal efficiency is obtained by the getCoor
function that takes two integers as arguments.
putCoor - assign coordinate values to a node. This function is
primarily aimed at preprocessor algorithms and should not be used
in moving grid problems (in the latter case one should apply the
move functions which are safer than putCoor, see the comments
above on deforming/moving grids).
getMinMaxCoord - returns the surrounding hypercube for the grid.
See class Grid.
getLocalCentroid - return the coordinates of the centroid of the
element in the local element coordinate system.
getGlobalCentroid - as getLocalCentroid, the coordinates are com
puted with respect to the global, physical coordinate system.
loc2glob - given a local node number in an element, the function
returns the corresponding global number of that node.
putLoc2Glob - used for setting the global node number of a local
node in an element.
setBoInd - set a particular boundary indicator at a node (ON
mode).
clearBoInd - turn off a particular boundary indicator at a node
(OFF mode).
putBoIndName - assign a name to a boundary indicator.
getBoIndName - get the name of a boundary indicator.
boSide - returns true if ALL nodes on a particular side in an
element are marked with a particular boundary indicator. The
function is often used for determining whether surface integrals
over an element are to be computed or not. However, the testing
inside the function is quite comprehensive. If the user is will
ing to store extra information, one can call
grid.additional_info.initSideIndicators (grid);
"GridFE::boSide will then use data structures in the addi
tional_info object (of type GridFEInfo) for a quick look up of
the side information.
boNode - returns true if a given node is marked with a particular
boundary indicator. There is an overloaded instance of this func
tion that returns true if the node is marked with any boundary
indicator.
boNodesInElm - returns true if any node in an element is marked
with a particular boundary indicator. There is an overloaded
instance of this function that returns true if any node in the
given element is marked with any boundary indicator.
BoInd - returns true if indicator ind is on in at least one node.
redefineBoInds - makes it possible to redefine the boundary indi
cators, that is, define a mapping between old and new boundary
indicators. Assume that we have three boundary indicators and
that we want to map these three indicators over to a new set of
four indicators. The new indicator 1 consists of nodes subjected
to old indicator 1 or 3, the new indicator 2 should be empty, the
new indicator three is identical to old indicator three, while
the new indicator four should consist of all the old indicators.
Let the names of the new indicators be "green, blue, red and
black. The mapping is defined by the following string which can
be sent as argument to redefine.
nbind=4, names= green blue red black 1=(1 3) 2=() 3=(3) 4=(1 2
3)
The names must be surrounded by blanks, only the two first equal
ity signs and the parenthesis are significant in the scan process
(thus the commands nbind, names, 1, 2 and 3 are never checked).
Old indicator numbers must be inside parenthesis and separated by
blanks. A equivalent string is
=4 =green blue red black (1 3)()(3)(1 2 3)
but less readable.
addBoIndNodes - makes it possible to mark new nodes with one of
the exisiting boundary indicators. New nodes are given as a
domain (hypercube), and the function will process all nodes that
are inside this domain. Here is an example of an input string to
the function.
"n=2 b1=[0,0]x[1,1] b4=[1,4]x[2,2]"
First the number of boundary indicators that are to be modified
(extended) must be given, it is 2 in this example. Then each
boundary indicator number of interest is preceeded by a "b" and
then its number, then an equal sign, and finally a tensor product
of intervals that defines the domain. All nodes inside this
domain ("inside" is determined by some tolerance, see Inter
val.inWIthTolerance, the tolerance equals the global variable
comparison_tolerance that is set in the initDiffpack function if
the command line argument +tolerance appears) are marked with the
boundary indicator. Here the point (0,1) will be marked with
indicator 1 and all nodes along the line , between and will be
marked with indicator 4. Sometimes one wants to mark just the
nodes inside the domain that are already marked by an indicator.
This is achieved by putting a capital B before the hypercubes:
"n=2 b1=B[0,0]x[1,1] b4=B[1,4]x[2,2]"
The B option is useful for curved boundaries since it enables
only the nodes on the curved boundary to be marked, and not
internal nodes as well.
The main application of this routine is when simple indication of
boundary indicators were made in the preprocessor, for example,
only the different sides of the domain were marked, and the user
wants to assign boundary conditions to single nodes or to parts
of a side. Often the redefineBoInds function is called first to
extend the number of indicators. The new indicators can then be
filled using addBoIndNodes. It can often be convenient to intro
duce an indicator that is not used in computations, but only for
marking the (whole) boundary. This makes it easier to use the
capital B option in front of the hypercubes (the boundary is then
already marked, and addBoIndNodes sets correct additional markers
for computational purposes). (In the input string, the first =,
the b and the = are the significant characters, strange errors
may occur if these characters are missing).
The final argument, bool add, can be given false or OFF as value.
In that case the indicators are turned off instead of on. Often
one must make two calls to addBoIndNodes, first turning some
indicators on and then turning other indicators off.
addMaterial - enables the user to give a hypercube and assign a
material number to all the elements that have their centroid
inside this hypercube. The syntax of the input string is:
"no=2 [0,2]x[-1,1]
where the first number after the equality sign is the material
number and the next string is the hypercube. In this particular
example, all elements that have the centroid inside will have
the material number set to 2. If the grid was made with only one
material (the default), all elements will have material number 1
initially.
In another version of addMaterial a polygon is used as input.
All elements with its centroid inside the polygon is set to the
given material number. The points defining the polygon is given
as a MatSimple matrix where the first and second columns define
the x and y coordinates of the polygon. The points in the polygon
should be given in a counter clockwise direction and the curve
has to be closed, i.e. the first and last element coincide. The
routine is mainly meant for 2D grids, but works also in 3D, the
polygon then defines a cylinder region. The addMaterial returns
the number of elements which are reset in the routine. inside
Polygon - tests if a point is inside the polygon. It returns 0
if the point is outside the polygon, 1 if inside and, 2 if the
point is on the border.
notEqual - determines whether two grids are equal. A safe proce
dure may be very time consuming (and is not implemented). The
parameter check_level reflects the amount of data that is com
pared. The value 1 implies a test on the address of the two grid
objects. Using such a criterion may result in a true return value
although the two grids are physically equal (they only correspond
to two copies). The default value is 2 for check_level which
means that the number of nodes and elements are compared. Future
extensions should detect the case where one grid is a deformation
of the other (differing nodal coordinates).
scanLattice - initializes the grid by reading a string with the
same syntax as GridLattice::scan. The grid will consist of
ElmTensorProd1 elements. The function is a convenient for gener
ating a finite element grid consisting of first order elements in
a box shaped domain, since there is no need to run a preproces
sor.
isLattice - returns true if the GridFE object is a regular lat
tice with uniform spacing and multi-linear elements (for example
produced by a PreproBox preprocessor). A true value from this
function makes it safe to call getLattice to get the correspond
ing GridLattice object that describeds the lattice. NOTE: If the
user has made the grid by the scanLattice function or the
PreproBox preprocessor, and later moved the coordinates using
putCoor, the grid is no longer a lattice, but isLattice will
return "true. Never use putCoor, try to use the safer move func
tions when transforming the grid coordinates. The move function
will ensure that isLattice will return false.
getLattice - returns access to a GridLattice object contained in
the GridFE object (if the finite element grid is regular with
uniform spacing). Call isLattice first to ensure that the grid
really is a lattice. If isLattice returns false, getLattice
returns actually a reference to a NULL pointer. What can you do
with the getLattice function? If you know that your finite ele
ment grid is actually a rectangular or box shaped grid with uni
form partition, you can operate on the finite element grid object
as if it were a finite difference type of grid. For example, the
spacing in the x-direction is easily available as grid.getLat
tice().Delta(1) if grid is type GridFE. Similarly, the number of
nodes in the y-direction is available as grid.getLattice().getDi
visions(2)+1.
isReallyLattice - performs a more thorough check that the grid is
really a lattice. Not implemented.
refineIfBox - makes a simple refinement of the grid if it was
created by a box preprocessor using multi-linear elements, or if
it was made by a call to scanLattice or operator=(const GridLat
tice&). The new grid is made by operator=(const GridLattice&).
This means that the element type will be, e.g., ElmB4n2D in 2D
even if the original type was ElmTensorProd1. The default bound
ary indicators should be intact (if they have been adjusted, a
warning is issued and one must remember to re-adjust the indica
tors). Note that this function gives an error message if the grid
is not a lattice (to avoid such messages, you can simply test on
"isLattice before you call refineIfBox). The argument to
refineIfBox reflects the refinement factor in each space
direction. For example, if the refinement factor is (2,2,2) in a
3D grid, the partition in each space direction is increased by a
factor of two (the mesh size is halved). For non-trivial refine
ments of finite element grids, we refer to class GridFEAdB.
setUniformMesh - sets an indicator for a uniform mesh to a true
value. By a uniform mesh we mean a mesh where all the elements
are of the same size and shape. In addition, the element type
(basis functions) should be identical for all elements. This is
often the case when a grid is computed by the PreproBox prepro
cessor. A uniform mesh where all the elements are of multi-linear
type is in our context denoted as a lattice grid. As an example,
a mesh consisting of biquadratic elements, all of the same size,
is uniform but the mesh is not a lattice. The setUniformMesh
function is usually used only in grid generating procedures.
setNonUniformMesh - sets the internal indicator for a uniform
mesh to a false value. The function also ensures that the grid is
not marked as a lattice.
hasUniformMesh - returns a true value if the grid is uniform. Our
definition of uniform is given in the documentation of the setU
niformMesh function.
move - displaces the grid according to a vector field (displace
ment) and a scale factor. Used for moving mesh problems. Can also
be used for plotting deformed mesh in, e.g., elasticity problems.
Then the scale factor should be determined from the size of the
grid - getMinMaxCoord - and the size of the maximum nodal dis
placement vector. Other overloaded versions of move enable trans
formations of the grid coordinates by means of a nodal coordinate
vector, another grid object and a vector field represented by
explicit functions. Always use move and not putCoor to transform
the coordinates of a grid after the grid is initially made.
startIterator - an iterator over all grid points is initialized
by this function.
nextPt - after having called startIterator one can call nextPt to
compute the "next" grid point. The grid point is given and
returned as an argument to the function. When nextPt returns
false there are no more points in the grid. The function is use
ful if it is desired to find the possible values of an explicit
function over the grid (loop through all grid points and compute
the maximum and minimum values).
getNoPoints - returns the number of nodes. It is fully equivalent
to "getNoNodes, but it is needed since the base class GridWithPts
specify getNoPoints as a virtual function.
isNode - returns a true value if the point given as argument is
really a node. If the grid is a lattice, a very efficient algo
rithm is run. In the general case, one calls GridWithPts::near
estPoint where all nodal points are compared to the given point.
If an equality (with a tolerance) is obtained, there is an exact
match and isNode returns the number of the actual node. If the
return value is 0, there is no exact match and the point is not a
node.
nearestPoint - a local variant of GridWithPts::nearestPoint. it
just calls that function. The reason to include nearestPoint as a
member function, is that GridWithPts defines the function to be
virtual, but provides no default version (a general static func
tion requiring a coordinate array for the grid points is offered
and it is this function that is called from GridFE:;nearest
Point).
isNodeInElm - returns a true value if the given node is a node in
the specified element.
boundingBox - calculates the max and min corrdinates of a hyper
cube, with axis parallell to the coordinate axis, that contains
the element in the interior.
findElmAndLocPt - given a global space point the routine finds
the element containing this point and the coordinates of the
point in the local element coordinate system. If the global point
is one of the nodes, the value of node is different from zero at
return (no element or local point are normally calculated in this
latter case). Sometimes one desires the element and local point
information even if the global point is a node. If the force_ele
ment_search variable is true, the element and local point will
always be found. If no element contains the global point, ele
ment=0 at return. The findElmAndLocPt function is typically used
in FieldFE::valuePt. There are two optional arguments for
returning the number of Newton iterations and the final "error"
of the inverse isoparametric mapping. The two arguments are
pointers and the data are returned only if the pointers are
nonzero at input. NULL is the default value of these pointer
arguments. One version of findElmAndLocPt calls a very efficient
algorithm in case of a lattice grid, and a standard search algo
rithm in the general case. The other version has a parameter Elm
SearchMethod that enables the programmer to choose among special
search methods. Since the search for elements and local points
can be very time consuming and also subject to problems, this
latter version of findElmAndLocPt is convenient for testing dif
ferent approaches. For normal application, the findElmAndLocPt
without the ElmSearchMethod parameter is the recommended version.
When the global point is not found in any element, the returned
values often contain DUMMY values - these will give rise to very
strange numerical results. Hence it should be easy to detect that
something wrong is happening. If one has trouble with
"FieldFE::valuePt, the problems are likely to be in findElmAnd
LocPt and it is recommended to re-run the case with the command
line option +verbose 6. This will result in messages from the
search algorithms. Problems will particularly be reported.
Brief documentation of the various search strategies are found in
the source file.
print - prints the grid.
scan - reads the grid from an input source.
newNodalNumbering - given an integer array old2new, this function
renumbers the nodes such that old node number i is assigned the
number old2new(i).
nodeCouplings, constraintCouplings - these functions determine
the coupling of nodes in the grid, that is, they constitute the
basis for computing the sparsity pattern (see the makeSparsity
Pattern routines in DegFreeFE.cpp).
calcBandwidth - computes the bandwidth associated with the mesh
(and one unknown at each node). We refer to the straightforward
source code for the precise definition of bandwidth in this con
text.
calcTolerance - a function that calculates the tolerance based on
the minimum element size in the grid divided by the argument
div_factor. The default value of div_factor is 100. The tolerance
is used when deciding if two nodes are equal. Note that the func
tion is static so it can be of general value when computing tol
erances in a grid.
The rest of the functions are usually only used by preprocessors.
We refer to the source code and the self-explanatory names for
information.
EXAMPLES
Normally, the computation of a finite element grid object of type
GridFE is carried out by a preprocessor. See reports on finite
element programming and preprocessors for examples. However, in
some occasions it is necessary to create a GridFE object based on
data from, e.g., a file. Next we will try to show how this can be
accomplished. That is, we will list all the functions that you
must call in order to initialize the GridFE object completely.
GridFE g;
Call the redim function of the GridFE object. This will allocate
the internal data structures. (Most cases will probably involve
nne=maxnne and one material.)
Is the grid actually a lattice? If so, one can rebind the lat
tice_data handle to a GridLattice object and then fill this
object with the proper information. Note that the lattice_data is
an option that is not required. If the lattice_data handle is
NULL (as it is by default) it just means that one cannot extract
lattice data like grid spacings in the different spatial direc
tions from a GridFE object. (Note that this point can be skipped
if you use unstructured or deformed grids.)
Uniform or non-uniform mesh: Call the setUniformMesh or the set
NonUniformMesh function. By default, the GridFE object assumes
that the grid is non uniform.
Set the element type (for the whole mesh, or for each element
individually) by calling setElmType. Example:
g.setElmType("ElmT4n3D"); // set type for all elements
Set the names of the boundary indicators by calling the put
BoIndName functions. Note that (i) the number of boundary indica
tors must correspond with the number given in the redim call, and
(ii) one can omit all boundary information in a GridFE object if
desired.
Run through all the nodes and call putCoor for each node to set
the nodal coordinates. Example:
g.putCoor(3.21, 2, 1); x-coord. at node 2 is 3.21
g.putCoor(3.22, 2, 2); y-coord. at node 2 is 3.22
Run through all the nodes and call setBoInd for the various
nodes that are subjected to boundary indicators. Example:
g.setBoInd (2, 5); indicator 5 at node 2 is on
Run through all the elements and compute the global node number
corresponding to the local nodes by calling the putLoc2Glob func
tion. In the same loop one should also set the subdomain number
of the element if there are more than one subdomains in the mesh.
Example:
g.putLoc2Glob (32, 3, 2); local node 2 in element 3 has global
number 32
g.putLoc2Glob (2, 2, 1);
g.putLoc2Glob (2, 5, 2);
g.putLoc2Glob (2, 4, 3);
Print a small grid to check the data. For larger grids one
should try to plot the grid and its boundaries.
SEE ALSO
class Grid, class ElmDef, class FieldFE
DEVELOPED BY
SINTEF Applied Mathematics, Oslo, Norway, and University of Oslo,
Dept. of Mathematics, Norway
AUTHOR
Hans Petter Langtangen, SINTEF/UiO. Some extensions by Klas
Samuelsson, UiO.