NAME
GridLattice - finite difference grid with constant partition
INCLUDE
include "GridLattice.h"
SYNTAX
class GridLattice : public GridWithPts
{
friend class GridFE;
friend class GridLatticeB;
friend class GridLatticeC;
friend class FieldLattice;
VecSimple(VecSimple(real)) coor; // coordinates in all space directions
int nsd; // number of space dimensions in grid (arbitrary)
Ptv(real) xmin; // min coordinate value in each spatial direction
Ptv(real) xmax; // max coordinate value in each spatial direction
Ptv(int) div; // number of grid intervals in each spatial direction
Ptv(real) scratch; // scratch array for internal computations
Ptv(int) base; // base value of each index
int npoints; // total number of lattice points
int ncells; // total number of lattice cells
// help variables for switching between single and multiple grid pt/box index
Ptv(int) multiple_idx;
Ptv(int) pidx_length; // for points, no of points in each space dir
Ptv(int) prod_pidx_length; // for points, products of pidx_length
Ptv(int) cidx_length; // for cells, no of cells in each space dir
Ptv(int) prod_cidx_length; // for cells, products of bidx_length
void init (); // compute idx_length, prod_idx_length, npoints etc
void initUniformCoordinates ();
void initUniformCoordinates (int i);
void initNonUniformCoordinates ();
// help variables for functions nextPt and startIterator :
Ptv(int) index_start_iterator, multidx, current_multidx, loop_length,
prod_loop_length;
int no_of_points;
int singleidx;
bool uniform; // true if gridspacing is uniform
// help variables for half indices (staggered grids)
// used in subclasses GridLatticeB, GridLatticeC
bool scalar_in_half_index; // true: scalar (p) in cell center
Ptv(real) index_base; // 0.5 if index in this dir is half,
// 0.0 if index is integer
Ptv(int) half_shift_in_min; // 1: half shift in this dir, 0: no half shift
Ptv(int) half_shift_in_max; // 1: half shift in this dir, 0: no half shift
protected:
int nbind; // number of boundary indicators
Handle(Indicators) bind; // boundary indicators and names
// help function for createIndexSet
IndexSet& findObstacles (const Ptv(int)& indicators,
const Ptv(int)& bmin,
const Ptv(int)& bmax);
public:
enum Gridtype { UNIFORM, NONUNIFORM };
enum Component { none=-1, p, u, v, w } comp_in_staggered_grid;
void setComponent (GridLattice::Component comp);
public:
GridLattice (int nsd = 1, enum Gridtype gtype = GridLattice:: UNIFORM);
~GridLattice ();
bool ok () const;
bool redim (int nsd);
bool notEqual (const GridLattice& grid) const;
bool isUniform () const { return uniform; }
void copy (const GridLattice& grid); // safer for derived classes
void operator = (const GridLattice& grid); // calls copy()
void initBoInds (); // sets up default boundary indicators
void initBoInds (int nbind); // only allocate space for 'nbind' bo inds
real& xMin (int i) { return xmin(i); }
real xMin (int i) const { return xmin(i); }
const Ptv(real)& xMin () const { return xmin; }
real& xMax (int i) { return xmax(i); }
real xMax (int i) const { return xmax(i); }
const Ptv(real)& xMax () const { return xmax; }
// grid spacing
real Delta (int dir) const; // spacing (for uniform grids only)
real Delta (int dir, int node) const;
const Ptv(real) Delta (const Ptv(int)& node) const;
real Delta (int dir, int node, int fict_boundary) const;
const Ptv(real) Delta (const Ptv(int)& node, int fict_boundary) const;
int getDivisions (int i) const // no of intervals in direction i
{ return div(i); }
const Ptv(int)& getDivisions () const
{ return div; }
Ptv(int) getNoGridPoints () const; // no of points in each direction
int getNoBoInds () const { return nbind; }
void fillCoords (int component, const VecSimple(real)& coords);
void fillSpacing (int component, const VecSimple(real)& spacing);
void setDivisions (int ndiv, int i); // not simple enough to allow int&
void setDivisions (const Ptv(int)& ndiv);
void setDivisions (const VecSimple(real)& coords, int i);
void setBase (int base_value) // base(i)=base_value, i=1,2,3,..
{ base.fill (base_value); }
int getBase (int i) const
{ return base(i); }
const Ptv(int)& getBase () const
{ return base; }
// index i runs from getBase(i) to getMaxI(i) in the grid:
int getMaxI (int i) const
{ return div(i) + base(i); }
Ptv(int) getMaxI () const;
const Ptv(int)& single2multiple (int single, bool points = true) const;
int multiple2single (const Ptv(int)& multiple,
bool points = true) const;
Ptv(real) getPt (const Ptv(int)& index) const; // general nD
real getPt (int dir, int index) const;
// note: getPt(single2multiple(i)), i is int, is possible!
Ptv(real) getPtCellCenter (const Ptv(int)& index) const;
real getPtCellCenter (int dir, int index) const;
void getMinMaxCellCoord (const Ptv(int)& index,
Ptv(real)& mincoord,
Ptv(real)& maxcoord) const;
// staggered grid
bool& scalarInHalfIndex() { return scalar_in_half_index; }
bool scalarInHalfIndex() const { return scalar_in_half_index; }
Ptv(int)& halfShiftInMax() { return half_shift_in_max; }
const Ptv(int)& halfShiftInMax() const { return half_shift_in_max; }
Ptv(int)& halfShiftInMin() { return half_shift_in_min; }
const Ptv(int)& halfShiftInMin() const { return half_shift_in_min; }
bool wholeIndex(real i) const;
bool halfIndex(real i) const;
virtual void setIndexBase();
virtual real& indexBase(int i) { return index_base(i); }
virtual real indexBase(int i) const { return index_base(i); }
virtual const Ptv(real)& indexBase () const { return index_base; }
// boundary indicator information:
String getBoIndName (int indno) const { return ( bind->getName(indno) ); }
bool boNode (int node, int indno) const { return bind->on(node,indno); }
bool boNode (int node) const; // is node marked with any indicator?
bool boInd (int ind); // any node with indicator no. ind?
void putBoIndName (const String& name, int indno)
{ bind->putName(name,indno); }
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); nbind = bind->nind; }
void addBoIndNodes (Is is, bool add = ON /*OFF: delete indicators*/);
// set functions mostly used by preprosessors
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); }
IndexSet& createIndexSetInnerNodes ();
IndexSet& createIndexSet (int ind);
IndexSet& createIndexSet (const VecSimple(int)& ind);
IndexSetCollection& createIndexSetCollection ();
bool insideDomain (const Ptv(real)& x, int fict_boundary=0) const;
bool insideDomain (const Ptv(int)& ind, int fict_boundary=0) const;
bool getCell (const Ptv(real)& x, Ptv(int)& index) const;
void findElmAndLocPt
(
const Ptv(real)& glob_pt,
int& element,
Ptv(real)& loc_pt,
int& node,
bool triangular_element = false
)const;
Ptv(int) nearestIndex
(const Ptv(real)& point, real& distance, bool& exact,
const int fict_boundary=0) const;
// virtual functions specified in the base classes:
virtual int nearestPoint
(const Ptv(real)& point, real& distance, bool& exact)
{ return multiple2single(nearestIndex(point,distance,exact)); }
virtual int getNoPoints () const;
virtual int getNoSpaceDim () const { return nsd; }
String print () const;
virtual void print (Os os) const;
virtual void scan (Is is);
virtual void scale ();
virtual void unscale ();
virtual void startIterator ();
void startIterator (const Ptv(int)& index1, const Ptv(int)& index2);
virtual bool nextPt (Ptv(real)& x);
bool nextPt (Ptv(int)& index);
bool nextPt (Ptv(real)& x, Ptv(int)& index);
virtual void getMinMaxCoord
(Ptv(real)& mincoord, Ptv(real)& maxcoord) const;
virtual void indexCheckOk(real i) const;
virtual void indexCheckOk(real i, real j) const;
virtual void indexCheckOk(real i, real j, real k ) const;
virtual void changeGridSize(){}
CLASS_INFO
VIRTUAL_CAST(GridLattice)
};
KEYWORDS
finite differences, grid, uniform, mesh
DESCRIPTION
GridLattice enables the programmer to define and manipulate a
rectangle/box shaped grid in any number of space dimensions. One
may specify UNIFORM or NONUNIFORM in the constructor, depending
on the required gridtype. A grid is only UNIFORM when it has uni
form partition in all space directions.
CONSTRUCTORS AND INITIALIZATION
The only constructor available has two arguments indicating the
size of internal arrays (ie. the number of space dimensions) and
the type of grid (UNIFORM or NONUNIFORM). The constructor can be
called with no arguments. The correct dimensions will anyway be
set by redim or by scan. The latter function is the usual tool to
initialize a GridLattice object. Here is an example:
GridLattice g;
g.scan ("d=2 [0,1]x[1,3] index:[0:20]x[0:56]");
All grid parameters must either be set by the scan function or
the programmer must use the redim function and then the functions
xMin, xMax and setDivisions (and perhaps setBase).
MEMBER FUNCTIONS
Most functions are self explanatory - see also the example below.
Note that there are internal scratch arrays (scratch and multi
ple_idx for example) that are used for returning references from
functions instead of copies. The single2multiple and getPt func
tions apply this technique. Be aware of this! There should be no
limitations on the use, except that new calls to these functions
will destroy the vector content that was returned in a previous
call.
The scan function reads an Is (string, stream, file) on the fol
lowing format for uniform partition
grid, d=3, [-1,1]x[0,1]x[0,1] indices=[0:10]x[1:8]x[-6:6]
This gives a 3D grid where -1 x 1, 0 y 1 and 0 z 1. The grid
points are numbered from 0 to 10 in the x-direction, hence there
are 11 points in the x-direction. In the y- and z-directions
there are 8 and 13 points respectively.
For a NONUNIFORM grid, the scan string is as follows:
grid, d=2 domain = [0,1]x[0,1] grid points = [0:5]x[1:10]
spacing=NONUNIFORM 0.0 0.2 0.4 0.6 0.8 1.0;
0.0 0.15 0.2 0.25 0.3 0.35 0.4 0.6 0.8 1.0;
The startIterator and nextPt functions enables the programmer to
run through all points in the grid or through only a subset of
the grid points. To initialize the iteration over points, call
startIterator. Without arguments, startIterator initializes the
iteration over all grid points. In an overloaded version indices
for the starting and stopping point in the grid can be given.
initBoInds - the version witch takes no argument sets up a set of
default boundary indicators along the edge of the grid, numbered
such that the first indicator is set on the top edge in the first
space-direction, the next at the top edge in the second space
direction and so on for all nsd, then the next indicator is set
on the lower edge in the first space direction, the next on the
lower edge of the second and so on. The version that takes an
integer argument does not set any indicators, it just allocates
space for the given number of boundary indicators, which must
later be set by the function addBoIndNodes.
xMin, xMax - gives the minimum and maximum coordinates of the
grid.
Delta - various versions of this function is implemented. The one
that takes one integer argument may only be used if the grid is
uniform and returns the spacing in the given direction. The ver
sion that takes two integer arguments returns the spacing
supplied node number. The version that takes a Ptv(int) argument,
is the same as the previous version, except it returns the
spacing in all space directions.
The two functions that also takes an integer specifying a fic
tional boundary is used if the FieldLattice object that uses this
grid has a ghost-boundary. This means that one is allowed to
request the spacing between nodes that either are inside the grid
or outside the grid, but inside the ghost-boundary. The spacings
returned for nodes in the ghost-boundary is the value of the cor
responding spacings on the grid, symmetrical to the edge of the
grid. See class FieldLattice for more info on ghost-boundaries.
getDivisions - returns the number of divisions ("cells") in a
given spatial direction.
getNoGridPoints - returns the number of grid points in each spa
tial direction.
getNoBoInds - returns the number of boundary indicators.
setDivisions - sets the number of grid points in one or all spa
tial direction depending on which version of the function is
called. This function then reinitializes the internal datamembers
that depend upon the number of divisions. The simplest version
of this function taking two integers can only be called for uni
form grids. The version that takes a Ptv(int) can be used for
both uniform and non-uniform grids, but if this function is
called for a non-uniform grid, it will initialize the coordinates
as if the grid was uniform. The default coordinates can then be
changed using fillCoords or "fillSpacings. The version that
takes a VecSimple(real) coordinate vector may be used for both
uniform and non-uniform grids, but it will set a uniform grid to
be non-uniform and use the supplied coordinates in the given
direction. An error message is issued if the supplied coordinates
do not coinside with the previously set xMin and xMax values.
fillCoords - fills the supplied coordinates into the datastruc
ture for the given space direction. The number of coordinates and
xMin and xMax values must coinside with the previously set val
ues.
fillSpacing - builds the coordinates in the given direction on
basis of the supplied spacings. Error messages are issued if the
number of spacings is incorrect or if the built coordinates do
not coinside with the previously set xMax value.
setBase - set the base of the indices counters (usually 0 or 1, 1
is default).
getBase - get the base of the indices counters.
getMaxI - returns the upper limit of a loop from the base to the
maximum number of grid points in a particular direction. For
example, in a 2D grid we may easily construct a loop over the
grid points:
int i0 = grid.getBase(1); int j0 = grid.getBase(2);
int in = grid.getMaxI(1); int jn = grid.getMaxI(2);
int i,j;
for (i = i0; i <= in; i++)
for (j = 1; j <= jn; j++)
// do something with point (i,j), e.g.
getPt - returns the coordinates of point. The point is specified
by its index. An overloaded version avoids the use of a Ptv(int)
for the index, instead one can give a direction and the index in
that direction.
getPtCellCenter - returns the coordinate of the center/centroid
of the cell corresponding to the given index. The numbering of
the cells starts at the base and ends with getMaxI minus one, in
each direction.
getCell - returns the d-dimensional index of the cell containing
the given spatial point.
nearestIndex - a la nearestPoint (see class GridWithPts), but
computes the index of the point that is nearest the specified
point. The index is an n-tuple (represented as usual by a
Ptv(int) object).
findElmAndLocPt - finds the cell and a local point in the cell
corresponding to a given spatial point. There are similar func
tions in class GridFE, but the present version is much more effi
cient since GridLattice is a grid with a very rigid geometric
structure. The cells in the GridLattice grid are treated as ten
sor product elements of first order (i.e., bilinear elements in
2D). There is also an optional argument for triangular elements
(only in 2D), the diagonal in the cell is then assumed to go from
the lower right to the upper left corner (this is compatible with
meshes generated by class PreproBox). The findElmAndLocPt func
tion first checks if the spatial point is a node, if it is the
corresponding element is chosen as the element that has the point
in the upper right corner. Hence the local element coordinates of
the point is typically (1,) (in 2D). If the spatial point is on
the lower boundary, x=xMin(i), the local point is of course
adjusted correspondingly. If the spatial point is not a node, one
computes the element and local coordinates from a formula. The
findElmAndLocPt function is automatically called from
GridFE::findElmAndLocPt if the finite element grid is actually a
lattice.
scalarInHalfIndex - true if the scalar (p) is located in the cen
triod of the grid cell in the staggered grid, false if the scalar
value is evaluated in the vertices of the cells. Where the vecoc
ity components are evaluated is determined by the evaluation
point for the scalar value and by the type of staggered grid (B-
grid or C-grid). See classes GridLatticeB and GridLatticeC for
more detailes.
halfShiftInMin - the half_shift_in_min variable is 0 by default,
i.e no half shift, a half shift is introduced by setting the
value 1 in the desired direction. The value must then be speci
fied in all directions by setting 0 in the directions that should
not have a half shift. This half shift removes the first half
grid cell in the given direction. This is convenient when we have
for instance the scalar value given at one end of the grid and a
velocity component given at the other end of the grid.
halfShiftInMax - same as halfShiftInMin but the half grid cell is
removed at the end of the grid in the given grid direction.
wholeIndex - true if the index i is an integer.
halfIndex - true if the index i is an integer plus one half.
setIndexBase - will usually not need to be called explicitly by
the user. Sets the index_base to the given values, depending on
the location of the various unknowns in the staggered grid. The
index_base is set such that if its value in the desired direction
is subtracted from the users supplied index, the resulting index
will always be an integer, allowing the correct array value to be
accessed. That is, if the staggered grid is such that a value is
evaluated in the center of the grid cell in a given direction,
the index_base is set to 0.5 in that same direction, otherwise
the index_base is set to 0.0.
indexBase - various functions returning the index base.
indexCheckOk - checks if the given indices are integers or inte
gers pluss one half, depending on the type of staggered grid
given.
changeGridSize - see classes GridLatticeB and GridLatticeC
insideDomain - checks if the given point (or the given indices)
is inside the grid. The optional argument fict_boundary should be
used when the corresponding FieldLattice object has a fictitious
boundary (ghost boundary) around the grid, and one wishes to find
out if the point/indices are inside the grid and the imaginary
boundary around the grid. See class FieldLattice for more info on
the effect of fictitious boundary (ghost boundary).
getBoIndName - get the name of the boundary indicator
boNode - the version with two arguments returns true if the given
node is marked with the given boundary indicator and false other
wise, and the other version of this function returns true if the
node is marked with at least one indicator.
boInd - returns true if at least one node in the grid is marked
with the given boundary indicator.
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.
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.
createIndexSetInnerNodes - creates an instance of the class
IndexSetBoxObstacles that contains all the nodes that does not
have any boundary indicators set.
createIndexSet - creates an instance of the class IndexSetBoxOb
stacles that contains all the nodes that has the specified/sup
plied boundary indicator(s) set, but none other than these set.
The function also fills the description string of the indexset
with information on which boundary indicators the indexset is
generated from.
createIndexSetCollection - creates IndexSetCollection that con
tains one indexset for each combination of boundary indicators,
plus one for the inner nodes. The IndexSetCollection is generated
by iterating through the grid with startIterator and nextPt. This
iteration is done such that i is the fastest varying index, then
j and finally k as the slowest varying index (in 3D).
Equivalently, in 2D, we have that for each index j, we run
through all the i-indices. (Here i, j and k are indices in direc
tion 1, 2 and 3 respectively.) The IndexSets organization are
determined by the order in which the boundary indicators are
found in the grid through this iteration. For each point in the
grid iteration, the boundary indicators set on this point are
examined. For each new combination found, a new indexset is added
to the collection, where this and future points that have this
particular combination are stored, hence the indexsets are stored
in the order the boundary indicator-combinations are found in the
grid. The indexset covering the inner nodes are always the last
indexset in the collection. Remember that one can always use the
getDescription function in each IndexSet for a description of
which boundary indicators that where the basis of the index set.
convertFromLattice - the functions are used for converting
IndexSetBoxObstacle objects to an IndexSetGen objects by the use
of the infomation from the supplied StencilCollection. (The func
tion is only available in the FDM module.)
See the base classes for documentation of most of the virtual
functions.
EXAMPLES
| // Define a finite difference grid in 2D on the rectangle
// [-1,1]x[0,10].
// Compute a sequence of grids where the cell length in each
space
// direction is reduced (to the half length f.ex.).
//
// In this example we will have
// xMin(1)=-1, xMin(2)=0, xMax(1)=1, xMax(2)=10
//
// Furthermore, the indices associated with grid starts at zero
(base(i)=0)
#include <GridLattice.h>
#include <IsOs.h>
int main (int nargs, const char** args)
{
initDiffpack (nargs, args); //not strictly necessary here, but
good habit
const int nsd = 2;
GridLattice grid(nsd);
grid.scan(cin);
int j;
// make 6 grids, reduction to half grid.Delta(i):
for (i = 1; i <= 6; i++)
{
for (j = 1; j <= nsd; j++)
grid.setDivisions (2*grid.getDivisions(j),j);
s_o << "Grid no. " << i;
s_o << " dx=" << grid.Delta(i) << " dy=" << grid.Delta(2)
<< ":0;
grid.print(s_o);
// what is the point with index (1,1) ?
s_o << " point (1,1) = ";
Ptv(int) index(2); index(1)=1; index(2)=1;
grid.getPt(index).print(s_o);
s_o << "0;
}
return 0;
}
The following input file should be prepared for the program
above:
------
Input grid, d=2 domain: [-1,1]x[0,10], indices [0:4,0:10]
------
SEE ALSO
class PartitionBox, class GeometryBox, class FieldLattice
DEVELOPED BY
SINTEF Applied Mathematics, Oslo, Norway, and University of Oslo,
Dept. of Mathematics, Norway
AUTHOR
Hans Petter Langtangen, SINTEF/UiO Extensions by Elizabeth Ack
lam, Numerical Objects AS