]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add two functions denoting measure and barycenter of a line and quad.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 6 Aug 1998 11:33:38 +0000 (11:33 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 6 Aug 1998 11:33:38 +0000 (11:33 +0000)
git-svn-id: https://svn.dealii.org/trunk@479 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/tria_accessor.h
deal.II/deal.II/source/grid/tria_accessor.cc

index f6c4365976f0e59d6a03f8de1b25bee977ccef37..3e5077d0d8ee5fee5b250b7dacdb21b5cceb8a22 100644 (file)
@@ -462,7 +462,23 @@ class LineAccessor :  public TriaAccessor<dim> {
                                      * midpoint of the line in real space.
                                      */
     Point<dim> center () const;
+
+                                    /**
+                                     * Return the barycenter of the line,
+                                     * which is the midpoint. The same
+                                     * applies as for the #center# function
+                                     * with regard to lines at the boundary.
+                                     */
+    Point<dim> barycenter () const;
     
+                                    /**
+                                     * Return the length of the line.
+                                     * The same
+                                     * applies as for the #center# function
+                                     * with regard to lines at the boundary.
+                                     */
+    double measure () const;
+
   private:
                                     /**
                                      *  Copy operator. This is normally
@@ -759,6 +775,27 @@ class QuadAccessor :  public TriaAccessor<dim> {
                                      */
     Point<dim> center () const;
 
+                                    /**
+                                     * Return the barycenter of the aud,
+                                     * which is the midpoint. The same
+                                     * applies as for the #center# function
+                                     * with regard to quads at the boundary.
+                                     */
+    Point<dim> barycenter () const;
+
+                                    /**
+                                     * Return the area of the quad. With
+                                     * area, we mean the area included by
+                                     * the straight lines connecting the
+                                     * four vertices, i.e. no information
+                                     * about the boundary is used; if
+                                     * you want other than bilinearly
+                                     * mapped unit quadrilaterals, ask the
+                                     * appropriate finite element class
+                                     * for the area of this quad.
+                                     */
+    double measure () const;
+
   private:
                                     /**
                                      *  Copy operator. This is normally
index a7b274bbaffd20b6c6fde37192959e6b05a9a543..953493289e0f3cd1ee645431d2adaebf8904e6fb 100644 (file)
@@ -307,6 +307,18 @@ Point<dim> LineAccessor<dim>::center () const {
 
 
 
+template <int dim>
+Point<dim> LineAccessor<dim>::barycenter () const {
+  return (vertex(1)+vertex(0))/2.;
+};
+
+
+
+template <int dim>
+double LineAccessor<dim>::measure () const {
+  return sqrt((vertex(1)-vertex(0)).square());
+};
+  
 
 
 
@@ -572,6 +584,131 @@ Point<dim> QuadAccessor<dim>::center () const {
 
 
 
+template <>
+Point<2> QuadAccessor<2>::barycenter () const {
+                                  // the evaluation of the formulae
+                                  // is a bit tricky when done dimension
+                                  // independant, so we write this function
+                                  // for 2D and 3D separately
+/*
+  Get the computation of the barycenter by this little Maple script. We
+  use the blinear mapping of the unit quad to the real quad. However,
+  every transformation mapping the unit faces to strait lines should
+  do.
+
+  Remember that the area of the quad is given by
+  |K| = \int_K 1 dx dy  = \int_{\hat K} |det J| d(xi) d(eta)
+  and that the barycenter is given by
+  \vec x_s = 1/|K| \int_K \vec x dx dy
+           = 1/|K| \int_{\hat K} \vec x(xi,eta) |det J| d(xi) d(eta)
+
+  # x and y are arrays holding the x- and y-values of the four vertices
+  # of this cell in real space. 
+  x := array(0..3);
+  y := array(0..3);
+  tphi[0] := (1-xi)*(1-eta):
+  tphi[1] := xi*(1-eta):
+  tphi[2] := xi*eta:
+  tphi[3] := (1-xi)*eta:
+  x_real := sum(x[s]*tphi[s], s=0..3):
+  y_real := sum(y[s]*tphi[s], s=0..3):
+  detJ := diff(x_real,xi)*diff(y_real,eta) - diff(x_real,eta)*diff(y_real,xi):
+
+  measure := simplify ( int ( int (detJ, xi=0..1), eta=0..1)):
+
+  xs := simplify (1/measure * int ( int (x_real * detJ, xi=0..1), eta=0..1)):
+  ys := simplify (1/measure * int ( int (y_real * detJ, xi=0..1), eta=0..1)):
+  readlib(C):
+
+  C(array(1..2, [xs, ys]), optimized);
+*/
+
+  const double x[4] = { vertex(0)(0),
+                       vertex(1)(0),
+                       vertex(2)(0),
+                       vertex(3)(0)  };
+  const double y[4] = { vertex(0)(1),
+                       vertex(1)(1),
+                       vertex(2)(1),
+                       vertex(3)(1)  };
+  const double t1 = x[0]*x[1];
+  const double t4 = x[2]*x[3];
+  const double t7 = x[0]*x[3];
+  const double t9 = x[1]*x[2];
+  const double t11 = x[0]*x[0];
+  const double t13 = x[1]*x[1];
+  const double t15 = x[2]*x[2];
+  const double t17 = x[3]*x[3];
+  const double t25 = -t1*y[0]+t1*y[1]-t4*y[2]+t4*y[3]+t7*y[0]-t9*y[1]
+                    +t11*y[1]-t13*y[0]+t15*y[3]-t17*y[2]-t11*y[3]
+                    +t13*y[2]-t15*y[1]+t17*y[0]+t9*y[2]-t7*y[3];
+  const double t35 = 1/(-x[0]*y[3]-x[1]*y[0]+x[2]*y[3]+x[3]*y[0]
+                       +x[0]*y[1]-x[3]*y[2]+x[1]*y[2]-x[2]*y[1]);
+  const double t38 = y[3]*x[3];
+  const double t40 = y[0]*x[0];
+  const double t42 = y[1]*x[1];
+  const double t44 = y[2]*x[2];
+  const double t50 = y[0]*y[0];
+  const double t52 = y[1]*y[1];
+  const double t54 = y[2]*y[2];
+  const double t56 = y[3]*y[3];
+  const double t62 = -t38*y[2]+t40*y[1]-t42*y[0]+t44*y[3]-t40*y[3]
+                    +t42*y[2]-t44*y[1]+t38*y[0]-x[1]*t50+x[0]*t52
+                    -x[3]*t54+x[2]*t56+x[3]*t50-x[2]*t52+x[1]*t54-x[0]*t56;
+  return Point<2> (t25*t35/3.0, t62*t35/3.0);
+};
+
+
+
+
+template <>
+double QuadAccessor<2>::measure () const {
+                                  // the evaluation of the formulae
+                                  // is a bit tricky when done dimension
+                                  // independant, so we write this function
+                                  // for 2D and 3D separately
+/*
+  Get the computation of the measure by this little Maple script. We
+  use the blinear mapping of the unit quad to the real quad. However,
+  every transformation mapping the unit faces to strait lines should
+  do.
+
+  Remember that the area of the quad is given by
+  \int_K 1 dx dy  = \int_{\hat K} |det J| d(xi) d(eta)
+
+  # x and y are arrays holding the x- and y-values of the four vertices
+  # of this cell in real space. 
+  x := array(0..3);
+  y := array(0..3);
+  tphi[0] := (1-xi)*(1-eta):
+  tphi[1] := xi*(1-eta):
+  tphi[2] := xi*eta:
+  tphi[3] := (1-xi)*eta:
+  x_real := sum(x[s]*tphi[s], s=0..3):
+  y_real := sum(y[s]*tphi[s], s=0..3):
+  detJ := diff(x_real,xi)*diff(y_real,eta) - diff(x_real,eta)*diff(y_real,xi):
+
+  measure := simplify ( int ( int (detJ, xi=0..1), eta=0..1)):
+  readlib(C):
+
+  C(measure, optimized);
+*/
+
+  const double x[4] = { vertex(0)(0),
+                       vertex(1)(0),
+                       vertex(2)(0),
+                       vertex(3)(0)  };
+  const double y[4] = { vertex(0)(1),
+                       vertex(1)(1),
+                       vertex(2)(1),
+                       vertex(3)(1)  };
+
+  return (-x[0]*y[3]/2.0-x[1]*y[0]/2.0+x[2]*y[3]/2.0+x[3]*y[0]/2.0+
+         x[0]*y[1]/2.0-x[3]*y[2]/2.0+x[1]*y[2]/2.0-x[2]*y[1]/2.0);
+};
+
+
+
 
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.