From: hartmann Date: Tue, 23 Dec 2008 13:13:17 +0000 (+0000) Subject: Bring patch deal.II_DLR@5587 to deal.II: Let the functions tranforming mother to... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=84736f2853b0e7a6505236b4e029104ef29fd451;p=dealii-svn.git Bring patch deal.II_DLR@5587 to deal.II: Let the functions tranforming mother to child coordinates (and vice versa) handle anisotropic children. git-svn-id: https://svn.dealii.org/trunk@18002 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/geometry_info.h b/deal.II/base/include/base/geometry_info.h index 5bb4175122..d03ed80bc4 100644 --- a/deal.II/base/include/base/geometry_info.h +++ b/deal.II/base/include/base/geometry_info.h @@ -1906,8 +1906,9 @@ struct GeometryInfo * depends on the number of the * child. */ - static Point cell_to_child_coordinates (const Point &p, - const unsigned int child_index); + static Point cell_to_child_coordinates (const Point &p, + const unsigned int child_index, + const RefinementCase refine_case=RefinementCase::isotropic_refinement); /** * The reverse function to the @@ -1917,8 +1918,9 @@ struct GeometryInfo * coordinate system of the * mother cell. */ - static Point child_to_cell_coordinates (const Point &p, - const unsigned int child_index); + static Point child_to_cell_coordinates (const Point &p, + const unsigned int child_index, + const RefinementCase refine_case=RefinementCase::isotropic_refinement); /** * Return true if the given point @@ -2357,16 +2359,271 @@ GeometryInfo::child_cell_from_point (const Point &) +template <> +inline +Point<1> +GeometryInfo<1>::cell_to_child_coordinates (const Point<1> &p, + const unsigned int child_index, + const RefinementCase<1> refine_case) + +{ + Assert (child_index < 2, + ExcIndexRange (child_index, 0, 2)); + Assert (refine_case==RefinementCase<1>::cut_x, + ExcInternalError()); + + return p*2.0-unit_cell_vertex(child_index); +} + + + +template <> +inline +Point<2> +GeometryInfo<2>::cell_to_child_coordinates (const Point<2> &p, + const unsigned int child_index, + const RefinementCase<2> refine_case) + +{ + Assert (child_index < GeometryInfo<2>::n_children(refine_case), + ExcIndexRange (child_index, 0, GeometryInfo<2>::n_children(refine_case))); + + Point<2> point=p; + switch (refine_case) + { + case RefinementCase<2>::cut_x: + point[0]*=2.0; + if (child_index==1) + point[0]-=1.0; + break; + case RefinementCase<2>::cut_y: + point[1]*=2.0; + if (child_index==1) + point[1]-=1.0; + break; + case RefinementCase<2>::cut_xy: + point*=2.0; + point-=unit_cell_vertex(child_index); + break; + default: + Assert(false, ExcInternalError()); + } + + return point; +} + + + +template <> +inline +Point<3> +GeometryInfo<3>::cell_to_child_coordinates (const Point<3> &p, + const unsigned int child_index, + const RefinementCase<3> refine_case) + +{ + Assert (child_index < GeometryInfo<3>::n_children(refine_case), + ExcIndexRange (child_index, 0, GeometryInfo<3>::n_children(refine_case))); + + Point<3> point=p; + // there might be a cleverer way to do this, but since this function is + // called in very few places for initialization purposes only, I don't + // care at the moment + switch (refine_case) + { + case RefinementCase<3>::cut_x: + point[0]*=2.0; + if (child_index==1) + point[0]-=1.0; + break; + case RefinementCase<3>::cut_y: + point[1]*=2.0; + if (child_index==1) + point[1]-=1.0; + break; + case RefinementCase<3>::cut_z: + point[2]*=2.0; + if (child_index==1) + point[2]-=1.0; + break; + case RefinementCase<3>::cut_xy: + point[0]*=2.0; + point[1]*=2.0; + if (child_index%2==1) + point[0]-=1.0; + if (child_index/2==1) + point[1]-=1.0; + break; + case RefinementCase<3>::cut_xz: + // careful, this is slightly + // different from xy and yz due to + // differnt internal numbering of + // children! + point[0]*=2.0; + point[2]*=2.0; + if (child_index/2==1) + point[0]-=1.0; + if (child_index%2==1) + point[2]-=1.0; + break; + case RefinementCase<3>::cut_yz: + point[1]*=2.0; + point[2]*=2.0; + if (child_index%2==1) + point[1]-=1.0; + if (child_index/2==1) + point[2]-=1.0; + break; + case RefinementCase<3>::cut_xyz: + point*=2.0; + point-=unit_cell_vertex(child_index); + break; + default: + Assert(false, ExcInternalError()); + } + + return point; +} + + + template inline Point -GeometryInfo::cell_to_child_coordinates (const Point &p, - const unsigned int child_index) +GeometryInfo::cell_to_child_coordinates (const Point &/*p*/, + const unsigned int /*child_index*/, + const RefinementCase /*refine_case*/) + { - Assert (child_index < GeometryInfo::max_children_per_cell, - ExcIndexRange (child_index, 0, GeometryInfo::max_children_per_cell)); + AssertThrow (false, ExcNotImplemented()); + return Point(); +} + + + +template <> +inline +Point<1> +GeometryInfo<1>::child_to_cell_coordinates (const Point<1> &p, + const unsigned int child_index, + const RefinementCase<1> refine_case) + +{ + Assert (child_index < 2, + ExcIndexRange (child_index, 0, 2)); + Assert (refine_case==RefinementCase<1>::cut_x, + ExcInternalError()); - return 2*p - unit_cell_vertex(child_index); + return (p+unit_cell_vertex(child_index))*0.5; +} + + + +template <> +inline +Point<3> +GeometryInfo<3>::child_to_cell_coordinates (const Point<3> &p, + const unsigned int child_index, + const RefinementCase<3> refine_case) + +{ + Assert (child_index < GeometryInfo<3>::n_children(refine_case), + ExcIndexRange (child_index, 0, GeometryInfo<3>::n_children(refine_case))); + + Point<3> point=p; + // there might be a cleverer way to do this, but since this function is + // called in very few places for initialization purposes only, I don't + // care at the moment + switch (refine_case) + { + case RefinementCase<3>::cut_x: + if (child_index==1) + point[0]+=1.0; + point[0]*=0.5; + break; + case RefinementCase<3>::cut_y: + if (child_index==1) + point[1]+=1.0; + point[1]*=0.5; + break; + case RefinementCase<3>::cut_z: + if (child_index==1) + point[2]+=1.0; + point[2]*=0.5; + break; + case RefinementCase<3>::cut_xy: + if (child_index%2==1) + point[0]+=1.0; + if (child_index/2==1) + point[1]+=1.0; + point[0]*=0.5; + point[1]*=0.5; + break; + case RefinementCase<3>::cut_xz: + // careful, this is slightly + // different from xy and yz due to + // differnt internal numbering of + // children! + if (child_index/2==1) + point[0]+=1.0; + if (child_index%2==1) + point[2]+=1.0; + point[0]*=0.5; + point[2]*=0.5; + break; + case RefinementCase<3>::cut_yz: + if (child_index%2==1) + point[1]+=1.0; + if (child_index/2==1) + point[2]+=1.0; + point[1]*=0.5; + point[2]*=0.5; + break; + case RefinementCase<3>::cut_xyz: + point+=unit_cell_vertex(child_index); + point*=0.5; + break; + default: + Assert(false, ExcInternalError()); + } + + return point; +} + + + +template <> +inline +Point<2> +GeometryInfo<2>::child_to_cell_coordinates (const Point<2> &p, + const unsigned int child_index, + const RefinementCase<2> refine_case) +{ + Assert (child_index < GeometryInfo<2>::n_children(refine_case), + ExcIndexRange (child_index, 0, GeometryInfo<2>::n_children(refine_case))); + + Point<2> point=p; + switch (refine_case) + { + case RefinementCase<2>::cut_x: + if (child_index==1) + point[0]+=1.0; + point[0]*=0.5; + break; + case RefinementCase<2>::cut_y: + if (child_index==1) + point[1]+=1.0; + point[1]*=0.5; + break; + case RefinementCase<2>::cut_xy: + point+=unit_cell_vertex(child_index); + point*=0.5; + break; + default: + Assert(false, ExcInternalError()); + } + + return point; } @@ -2374,16 +2631,16 @@ GeometryInfo::cell_to_child_coordinates (const Point &p, template inline Point -GeometryInfo::child_to_cell_coordinates (const Point &p, - const unsigned int child_index) +GeometryInfo::child_to_cell_coordinates (const Point &/*p*/, + const unsigned int /*child_index*/, + const RefinementCase /*refine_case*/) { - Assert (child_index < GeometryInfo::max_children_per_cell, - ExcIndexRange (child_index, 0, GeometryInfo::max_children_per_cell)); - - return (p + unit_cell_vertex(child_index))/2; + AssertThrow (false, ExcNotImplemented()); + return Point(); } + template <> inline bool