From: Wolfgang Bangerth Date: Thu, 6 Aug 1998 11:33:38 +0000 (+0000) Subject: Add two functions denoting measure and barycenter of a line and quad. X-Git-Tag: v8.0.0~22775 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7e9574b7fd2022ce14aa2693a66ad14d69520c25;p=dealii.git Add two functions denoting measure and barycenter of a line and quad. git-svn-id: https://svn.dealii.org/trunk@479 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/grid/tria_accessor.h b/deal.II/deal.II/include/grid/tria_accessor.h index f6c4365976..3e5077d0d8 100644 --- a/deal.II/deal.II/include/grid/tria_accessor.h +++ b/deal.II/deal.II/include/grid/tria_accessor.h @@ -462,7 +462,23 @@ class LineAccessor : public TriaAccessor { * midpoint of the line in real space. */ Point 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 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 { */ Point 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 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 diff --git a/deal.II/deal.II/source/grid/tria_accessor.cc b/deal.II/deal.II/source/grid/tria_accessor.cc index a7b274bbaf..953493289e 100644 --- a/deal.II/deal.II/source/grid/tria_accessor.cc +++ b/deal.II/deal.II/source/grid/tria_accessor.cc @@ -307,6 +307,18 @@ Point LineAccessor::center () const { +template +Point LineAccessor::barycenter () const { + return (vertex(1)+vertex(0))/2.; +}; + + + +template +double LineAccessor::measure () const { + return sqrt((vertex(1)-vertex(0)).square()); +}; + @@ -572,6 +584,131 @@ Point QuadAccessor::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); +}; + + +