Toast++  1.0.2 (r.539)
Forward and inverse modelling in optical tomography
tri3.h
1 // -*-C++-*-
2 // ==========================================================================
3 // Module libfe
4 // File tri3.h
5 // Declaration of unstructured 2-D element class Triangle3
6 // (3-noded triangle, 1st order shape functions)
7 //
8 //
9 // Node arrangement: N2
10 // ----------------- + The local element has node coordinates
11 // |\ N0 = (0,0)
12 // | \ N1 = (1,0)
13 // | \ N2 = (0,1)
14 // | \
15 // | \ Sides: side 0 contains N0,N1 (y=0)
16 // | \ side 1 contains N1,N2 (x+y=1)
17 // | \ side 2 contains N2,N0 (x=0)
18 // N0+-------+N1
19 //
20 // Inheritance:
21 // ------------
22 // Element
23 // ---> Element_Unstructured
24 // ---> Element_Unstructured_2D
25 // ---> Triangle3
26 // ==========================================================================
27 
28 #ifndef __TRI3_H
29 #define __TRI3_H
30 
31 #include "toastdef.h"
32 
33 class Surface;
34 class Mesh;
35 
42 class FELIB Triangle3: public Element_Unstructured_2D {
43 public:
44 
48  Triangle3();
49 
54  Triangle3 (const Triangle3& el);
55 
59  ~Triangle3();
60 
64  Element *Copy();
65 
74  void Initialise (const NodeList& nlist);
75 
76  BYTE Type () const { return ELID_TRI3; }
77  // returns element type id
78 
79  BYTE VtkType() const { return 5; }
80 
81  unsigned long GetCaps () const { return ELCAPS_SUBSAMPLING; }
82  // Return element capability flags
83 
84  int nNode() const { return 3; }
85  // returns number of nodes per element
86 
87  int nSide () const { return 3; };
88  // number of sides per element
89 
90  int nSideNode (int /*side*/) const { return 2; };
91  // returns number of nodes attached to side 'side'
92 
93  int SideNode (int side, int node) const;
94  // returns number of 'node'-th node of side 'side'
95 
96  double SideSize (int sd, const NodeList &nlist) const;
97  // returns length of side 'sd'
98 
105  double DetJ (const Point &loc, const NodeList *nlist = 0) const
106  { return Size()*2.0; }
107 
108  Point Local (const NodeList& nlist, const Point& glob) const;
109  // returns the local coordinate corresponding to global coordinate 'glob'
110 
111  Point NodeLocal (int node) const;
112  // returns local coordinates of node 'node'
113 
114  Point SurfToLocal (int side, const Point &p) const;
115  // Map 1-D surface point p (0-1) along side 'side' into 2-D local coords
116 
117  RVector DirectionCosine (int side, RDenseMatrix &jacin);
118  // returns a vector of direction cosines in global coordinates of the
119  // normal to side 'side'.
120 
121  const RVector &LNormal (int side) const;
122 
123  bool LContains (const Point& loc, bool pad = true) const;
124  // returns TRUE if point 'loc' (in local coordinates of the element) is
125  // within the standard triangle
126 
127  bool GContains (const Point& glob, const NodeList& nlist) const;
128  // returns TRUE if point 'glob' (in global coordinates) is within the
129  // element. nlist is the mesh node coordinate list
130 
131  RVector LocalShapeF (const Point &loc) const;
132  // vector of shape functions u_i(loc) at point 'loc' given in local
133  // element coordinates
134 
135  RDenseMatrix LocalShapeD (const Point &loc) const;
136  // 2x3 matrix of shape function derivatives (d u_j)/(d x_i) at point 'loc'
137  // given in local element coordinates
138 
139  RVector GlobalShapeF (const NodeList& nlist, const Point& glob) const;
140  // returns the shape functions for each node at a point, given in global
141  // coordinates
142 
143  RDenseMatrix GlobalShapeD (const NodeList& nlist, const Point& glob) const;
144  // returns derivatives of shape functions in global coordinates at point
145  // `glob', given in global coordinates
146 
147  double IntF (int i) const;
148 
149  RSymMatrix IntFF () const;
150  // Return integral over element of product of shape functions:
151  // FF = Int_el { F_i(r) F_j(r) } dr
152 
153  int QuadRule (int order, const double **wght, const Point **absc) const;
154 
155  double IntFF (int i, int j) const;
156  // Return a single element of IntFF
157 
158  double IntFFF (int i, int j, int k) const;
159  // returns a single element of integral over element of product of three
160  // shape functions:
161  // IntFFF = Int_el { F_i(r) * F_j(r) * F_k(r) } dr
162 
163 #ifdef UNDEF
164  // MS 29.6.99 Removed because incompatible with quadratic shape functions
165 
166  void IntFFF (double &iii, double &iij, double &ijk) const;
167  // returns the values of the FFF tensor for:
168  // all indices equal (iii)
169  // two indices equal (iij)
170  // all indices different (ijk)
171 #endif
172 
173  RSymMatrix IntPFF (const RVector& P) const;
174  // Returns integral over element of product of two shape functions and a
175  // function P defined on nodes:
176  // PFF = Int_el { P(r) F_i(r) F_J(r) } dr
177  // Vector P contains the complete nodal solution for the mesh
178 
179  double IntPFF (int i, int j, const RVector& P) const;
180  // Returns a single element of IntPFF
181 
182  RVector IntFD (int i, int j) const;
183  // Returns single (vector) element of FD integral over element:
184  // IntFD = Int_el { F_i(r) D_j(r) } dr
185 
186  double IntFDD (int i, int j, int k) const
187  { return intdd.Get(j,k) * 0.3333333333; }
188  // returns a single element of integral over element:
189  // IntFDD = Int_el { F_i(r) * D_j(r) * D_k(r) } dr
190 
191  RSymMatrix IntPDD (const RVector& P) const
192  { return intdd * ((P[Node[0]]+P[Node[1]]+P[Node[2]])/3.0); }
193  // Returns integral over element of product of two shape derivatives and a
194  // function P defined in nodal basis:
195  // PDD = Int_el { P(r) D_i(r) D_j(r) } dr
196  // Vector P contains the complete nodal solution for the mesh
197 
198  double IntPDD (int i, int j, const RVector &P) const
199  { return intdd.Get(i,j) * ((P[Node[0]]+P[Node[1]]+P[Node[2]])/3.0); }
200  // Returns a single element of IntPDD
201 
202  RSymMatrix BndIntPFF (const RVector &P) const
203  { return intbff * ((P[Node[0]]+P[Node[1]]+P[Node[2]])/3.0); }
204  // This is a hack! Really P would have to be taken into the integral,
205  // rather than just averaged.
206 
207  double BndIntPFF (int i, int j, const RVector &P) const
208  { return intbff.Get(i,j) * ((P[Node[0]]+P[Node[1]]+P[Node[2]])/3.0); }
209  // Returns a single element of BndIntPFF
210 
211  double IntFd (int i, int j, int k) const;
212  // Returns a single element of Int [u_i du_j/dx_k] dr
213 
214  double IntPd (const RVector &P, int j, int k) const;
215  // Returns Int [p(r) du_j/dx_k] dr where p(r) is defined by its nodal
216  // basis coefficients P.
217 
218  double IntFdd (int i, int j, int k, int l, int m) const;
219  // Int u_i du_j/dx_l du_k/dx_m dr
220 
221  double IntPdd (const RVector &p, int j, int k, int l, int m) const;
222  // Int f(r) du_j/dx_l du_k/dx_m dr
223  // where f(r) is given as a nodal vector
224 
225  RSymMatrix Intdd() const;
226  // returns matrix of mixed derivatives
227 
228  double BndIntFSide (int i, int sd);
229 
230  double BndIntFFSide (int i, int j, int sd);
231  // Int [u_i u_j] dr along side sd
232 
233  RVector BndIntFX (int side, double (*func)(const Point&),
234  const NodeList &nlist) const;
235  // See element.h for description
236 
237  RVector BndIntFCos (int side, const Surface *surf, const RVector &cntcos,
238  double sigma, double sup, const NodeList &nlist) const;
239  // See element.h for description
240 
241  RVector BndIntFCos (int side, const RVector &cntcos, double a,
242  const NodeList &nlist) const;
243  // See element.h for description
244 
245  RVector BndIntFDelta (int side, const Surface *surf, const RVector &pos,
246  const NodeList &nlist) const;
247  // See element.h for description
248 
249  int GetLocalSubsampleAbsc (const Point *&absc) const;
250  // returns abscissae in local coords for numerical integration by uniform
251  // subsampling. Return value is the number of samples returned in 'absc'
252 
253  int GetBndSubsampleAbsc (int side, const Point *&absc) const;
254  // abscissae for numerical integration over boundary side in local
255  // coordinates of boundary element (dim-1). Use SurfToLocal to convert
256  // into local element coordinates
257 
258  int GlobalIntersection (const NodeList &nlist, const Point &p1,
259  const Point &p2, Point **list);
260  // Same as 'Intersection' but p1 and p2 given in global coords
261  // The return list however is still in local coords
262 
263  int Intersection (const Point &p1, const Point &p2, Point** pi);
264  // creates a list of points where the line defined by p1 and p2 intersects
265  // the element (in local coordinates) or starts/ends within the element.
266  // Returns the length of the list
267 
268  void SplitSide (Mesh *mesh, int side, int newnode,
269  Element *nbr1, Element *nbr2, Element *el1, Element *el2);
270  // Subdivide the triangle such that 'side' is split into two parts.
271  // If newnode >= 0 and nbr1, nbr2 != NULL, then the side neighbour has
272  // already been split. In that case, newnode is the node index of the
273  // midpoint node, and nbr1, nbr2 are the new neighbours for the two
274  // side segments.
275  // Otherwise, the SplitSide method is also called for the neighbour
276  // SplitSide either performs a bisection if the current subdivision
277  // level is even, otherwise it performs a merge and resplit operation.
278  // In the latter case, it also calls SplitSide for any additionally
279  // subdivided sides.
280  // On return, el1 and el2 are pointers to the two elements that now
281  // face the split side (i.e. the neighbours of nbr1 and nbr2, if those
282  // were provided)
283 
284 
285  void Bisect (Mesh *mesh, int side, int newnode,
286  Element *nbr1, Element *nbr2, Element *el1, Element *el2);
287 
288  void MergeAndResplit (Mesh *mesh, int side, int newnode,
289  Element *nbr1, Element *nbr2, Element *el1, Element *el2);
290 
291 protected:
292 
293  double ComputeSize (const NodeList &nlist) const;
294  // area of triangle in global coordinates
295 
296  RSymMatrix ComputeIntDD (const NodeList &nlist) const;
297  // Returns integral over element of product of shape derivatives:
298  // DD = Int_el { D_i(r) D_j(r) } dr
299 
300  void ComputeIntFD (const NodeList &nlist);
301  // Calculates the FD integrals over the element:
302  // IntFD = Int_el { F_i(r) D_j(r) } dr
303 
304  RSymMatrix ComputeBndIntFF (const NodeList& nlist) const;
305  // Returns line integral of product of two shape functions along sides of
306  // the element which belong to the mesh boundary
307  // If the element does not contain any boundary sides then the returned
308  // matrix is empty.
309 
310  //private:
311 
312  double *intfd_0, *intfd_1;
313  // IntFD matrix storage
314 
315  double a0, b0, c0, a1, b1, c1, a2, b2, c2;
316  // triangle geometry parameters
317 
318 #ifdef TRI3_STORE_INTFF
319  RSymMatrix intff;
320  // stores integral over element of product of shape functions
321  // Int_el { F_i(r) F_j(r) } dr
322  // set by Initialise
323 #endif
324 
325 };
326 
327 #endif
Definition: ndlist.h:14
Templated vector class.
Definition: vector.h:39
virtual bool GContains(const Point &glob, const NodeList &nlist) const
Checks if a global point coordinate is inside the element.
virtual RSymMatrix IntFF() const =0
Integrals of all products of two shape functions over the element.
double DetJ(const Point &loc, const NodeList *nlist=0) const
Returns determinant of Jacobian.
Definition: tri3.h:105
double Size() const
Returns the element size.
Definition: element.h:1128
Definition: node.h:39
virtual double IntFdd(int i, int j, int k, int l, int m) const
Integral of the product of a shape function and two partial shape function derivatives over the eleme...
Definition: element.h:850
Definition: surface.h:21
virtual double IntPd(const RVector &P, int j, int k) const
Integral of the product of a nodal function and a partial shape function derivative over the element...
Definition: element.h:823
BYTE VtkType() const
Returns the VTK element type identifier, or 0 if the element doesn't have a VTK representation.
Definition: tri3.h:79
#define ELID_TRI3
3-noded triangle
Definition: element.h:45
Finite-element mesh management.
Definition: mesh.h:145
virtual RVector DirectionCosine(int side, RDenseMatrix &jacin)=0
Returns the direction cosines of a side normal.
virtual RVector GlobalShapeF(const NodeList &nlist, const Point &glob) const
Returns the values of the shape functions at a global point.
Definition: element.h:490
int nNode() const
Returns the number of nodes associated with the element.
Definition: tri3.h:84
virtual double IntFFF(int i, int j, int k) const =0
Integral of a product of three shape functions over the element.
Definition: point.h:18
virtual double BndIntFFSide(int i, int j, int sd)=0
Surface integral of a product of two shape functions over one of the sides of the element...
virtual RDenseMatrix LocalShapeD(const Point &loc) const =0
Returns the values of the shape function derivatives at a local point.
double IntPDD(int i, int j, const RVector &P) const
Integrals of a product of a nodal function and two shape function derivatives over the element...
Definition: tri3.h:198
double BndIntPFF(int i, int j, const RVector &P) const
Surface integrals of a product of a nodal function and two shape functions over all boundary sides of...
Definition: tri3.h:207
virtual double IntFd(int i, int j, int k) const
Integral of the product of a shape function and a partial shape function derivative over the element...
Definition: element.h:810
virtual const RVector & LNormal(int side) const =0
Returns a side normal in local coordinates.
virtual RSymMatrix IntPFF(const RVector &P) const =0
Integrals of all products of two shape functions and a nodal function over the element.
virtual double BndIntFSide(int i, int sd)
Surface integral of a shape function over one of the sides of the element.
Definition: element.h:707
virtual void Initialise(const NodeList &nlist)
Element initialisation.
virtual RDenseMatrix GlobalShapeD(const NodeList &nlist, const Point &glob) const
Returns the values of the shape function derivatives at a global point.
Definition: element.h:502
A 3-noded (linear) 2-dimensional triangle element.
Definition: tri3.h:42
double IntFDD(int i, int j, int k) const
Integral of a product of a shape function and two shape function derivatives over the element...
Definition: tri3.h:186
unsigned long GetCaps() const
Returns element capability flags.
Definition: tri3.h:81
virtual RVector IntFD(int i, int j) const
Integral of a product of a shape function and a shape function derivative over the element...
Definition: element.h:637
Base class for finite element types.
Definition: element.h:84
BYTE Type() const
Returns an element type identifier.
Definition: tri3.h:76
Base class for all 2-D unstructured element types.
Definition: element.h:1234
virtual double IntF(int i) const =0
Integral of a shape function over the element.
int nSide() const
Returns the number of element sides.
Definition: tri3.h:87
virtual double SideSize(int side, const NodeList &nlist) const
Returns the size of an element side.
Definition: element.h:397
virtual RVector LocalShapeF(const Point &loc) const =0
Returns the values of the shape functions at a local point.
RSymMatrix IntPDD(const RVector &P) const
All integrals of products of a nodal function and two shape function derivatives over the element...
Definition: tri3.h:191
virtual Point SurfToLocal(int side, const Point &p) const
Maps a point from surface coordinates to local element coordinates.
Definition: element.h:302
virtual RSymMatrix Intdd() const
Integral of the product of two partial shape function derivatives over the element.
Definition: element.h:836
virtual Element * Copy()=0
Create a copy of the element and return a pointer to it.
Dense matrix class.
Definition: crmatrix.h:38
virtual Point Local(const NodeList &nlist, const Point &glob) const =0
Maps a point from global to local element coordinates.
RSymMatrix BndIntPFF(const RVector &P) const
Surface integrals of all products of a nodal function and two shape functions over all boundary sides...
Definition: tri3.h:202
virtual Point NodeLocal(int node) const =0
Returns the local coordinates of an element node.
#define ELCAPS_SUBSAMPLING
element implements IntFFF and IntFDD by subsampling
Definition: element.h:61
virtual int SideNode(int side, int node) const =0
Returns relative node index for a side vertex.
virtual int QuadRule(int order, const double **wght, const Point **absc) const
Returns the weights and abscissae of quadrature rules over the element.
Definition: element.h:521
virtual bool LContains(const Point &loc, bool pad=true) const =0
Checks if a local point coordinate is inside the element.
int nSideNode(int) const
Returns the number of vertices associated with a side.
Definition: tri3.h:90