From e192f4ef5dd241e13567903d1678377929178a92 Mon Sep 17 00:00:00 2001 From: wolf Date: Sat, 20 Feb 1999 14:48:30 +0000 Subject: [PATCH] Move to the new way of describing boundaries. git-svn-id: https://svn.dealii.org/trunk@862 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/Attic/examples/poisson/problem.cc | 128 ++++++++++++++---- tests/big-tests/poisson/problem.cc | 128 ++++++++++++++---- 2 files changed, 202 insertions(+), 54 deletions(-) diff --git a/deal.II/deal.II/Attic/examples/poisson/problem.cc b/deal.II/deal.II/Attic/examples/poisson/problem.cc index c683811385..20ffc09441 100644 --- a/deal.II/deal.II/Attic/examples/poisson/problem.cc +++ b/deal.II/deal.II/Attic/examples/poisson/problem.cc @@ -16,7 +16,11 @@ class BoundaryValuesSine : public Function { * at the given point. */ virtual double operator () (const Point &p) const { - return cos(2*3.1415926536*p(0))*cos(2*3.1415926536*p(1)); + double x = 1; + + for (unsigned int i=0; i { Assert (values.size() == points.size(), ExcVectorHasWrongSize(values.size(), points.size())); for (unsigned int i=0; i::operator() (points[i]); }; }; @@ -100,25 +103,96 @@ template class CurvedLine : public StraightBoundary { public: - virtual Point in_between (const PointArray &neighbors) const { - Point middle = StraightBoundary::in_between(neighbors); - double x=middle(0), - y=middle(1); - - if (y + get_new_point_on_line (const typename Triangulation::line_iterator &line) const; - return middle; - }; + virtual Point + get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; +}; + + + +template +Point +CurvedLine::get_new_point_on_line (const typename Triangulation::line_iterator &line) const +{ + Point middle = StraightBoundary::get_new_point_on_line (line); + + // if the line is at the top of bottom + // face: do a special treatment on + // this line. Note that if the + // z-value of the midpoint is either + // 0 or 1, then the z-values of all + // vertices of the line is like that + if (dim>=3) + if (((middle(2) == 0) || (middle(2) == 1)) + // find out, if the line is in the + // interior of the top or bottom face + // of the domain, or at the edge. + // lines at the edge need to undergo + // the usual treatment, while for + // interior lines taking the midpoint + // is sufficient + // + // note: the trick with the boundary + // id was invented after the above was + // written, so we are not very strict + // here with using these flags + && (line->boundary_indicator() == 1)) + return middle; + + + double x=middle(0), + y=middle(1); + + if (y +Point +CurvedLine::get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const +{ + Point middle = StraightBoundary::get_new_point_on_quad (quad); + + // if the face is at the top of bottom + // face: do not move the midpoint in + // x/y direction. Note that if the + // z-value of the midpoint is either + // 0 or 1, then the z-values of all + // vertices of the quad is like that + if ((middle(2) == 0) || (middle(2) == 1)) + return middle; + + double x=middle(0), + y=middle(1); + + if (y::make_zoom_in_grid () { tria->execute_refinement (); Triangulation::active_cell_iterator cell; - for (int i=0; i<17; ++i) + for (int i=0; i<(dim==3 ? 5 : 17); ++i) { // refine the presently - // second last cell 17 + // second last cell several // times cell = tria->last_active(tria->n_levels()-1); --cell; @@ -362,7 +436,7 @@ void PoissonProblem::make_random_grid () { tria->refine_global (1); Triangulation::active_cell_iterator cell, endc; - for (int i=0; i<12; ++i) + for (int i=0; i<(dim==3 ? 7 : 12); ++i) { int n_levels = tria->n_levels(); cell = tria->begin_active(); @@ -455,7 +529,7 @@ void PoissonProblem::run (ParameterHandler &prm) { if (!set_boundary_values (prm)) return; - FELinear fe; + FEQuadraticSub fe; PoissonEquation equation (*rhs); QGauss2 quadrature; @@ -481,7 +555,7 @@ void PoissonProblem::run (ParameterHandler &prm) { string o_filename = prm.get ("Output file"); ofstream gnuplot(o_filename.c_str()); fill_data (out); - out.write_gnuplot (gnuplot); + out.write_ucd (gnuplot); gnuplot.close (); // release the lock of the DoF object to @@ -495,4 +569,4 @@ void PoissonProblem::run (ParameterHandler &prm) { -template class PoissonProblem<2>; +template class PoissonProblem<3>; diff --git a/tests/big-tests/poisson/problem.cc b/tests/big-tests/poisson/problem.cc index c683811385..20ffc09441 100644 --- a/tests/big-tests/poisson/problem.cc +++ b/tests/big-tests/poisson/problem.cc @@ -16,7 +16,11 @@ class BoundaryValuesSine : public Function { * at the given point. */ virtual double operator () (const Point &p) const { - return cos(2*3.1415926536*p(0))*cos(2*3.1415926536*p(1)); + double x = 1; + + for (unsigned int i=0; i { Assert (values.size() == points.size(), ExcVectorHasWrongSize(values.size(), points.size())); for (unsigned int i=0; i::operator() (points[i]); }; }; @@ -100,25 +103,96 @@ template class CurvedLine : public StraightBoundary { public: - virtual Point in_between (const PointArray &neighbors) const { - Point middle = StraightBoundary::in_between(neighbors); - double x=middle(0), - y=middle(1); - - if (y + get_new_point_on_line (const typename Triangulation::line_iterator &line) const; - return middle; - }; + virtual Point + get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; +}; + + + +template +Point +CurvedLine::get_new_point_on_line (const typename Triangulation::line_iterator &line) const +{ + Point middle = StraightBoundary::get_new_point_on_line (line); + + // if the line is at the top of bottom + // face: do a special treatment on + // this line. Note that if the + // z-value of the midpoint is either + // 0 or 1, then the z-values of all + // vertices of the line is like that + if (dim>=3) + if (((middle(2) == 0) || (middle(2) == 1)) + // find out, if the line is in the + // interior of the top or bottom face + // of the domain, or at the edge. + // lines at the edge need to undergo + // the usual treatment, while for + // interior lines taking the midpoint + // is sufficient + // + // note: the trick with the boundary + // id was invented after the above was + // written, so we are not very strict + // here with using these flags + && (line->boundary_indicator() == 1)) + return middle; + + + double x=middle(0), + y=middle(1); + + if (y +Point +CurvedLine::get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const +{ + Point middle = StraightBoundary::get_new_point_on_quad (quad); + + // if the face is at the top of bottom + // face: do not move the midpoint in + // x/y direction. Note that if the + // z-value of the midpoint is either + // 0 or 1, then the z-values of all + // vertices of the quad is like that + if ((middle(2) == 0) || (middle(2) == 1)) + return middle; + + double x=middle(0), + y=middle(1); + + if (y::make_zoom_in_grid () { tria->execute_refinement (); Triangulation::active_cell_iterator cell; - for (int i=0; i<17; ++i) + for (int i=0; i<(dim==3 ? 5 : 17); ++i) { // refine the presently - // second last cell 17 + // second last cell several // times cell = tria->last_active(tria->n_levels()-1); --cell; @@ -362,7 +436,7 @@ void PoissonProblem::make_random_grid () { tria->refine_global (1); Triangulation::active_cell_iterator cell, endc; - for (int i=0; i<12; ++i) + for (int i=0; i<(dim==3 ? 7 : 12); ++i) { int n_levels = tria->n_levels(); cell = tria->begin_active(); @@ -455,7 +529,7 @@ void PoissonProblem::run (ParameterHandler &prm) { if (!set_boundary_values (prm)) return; - FELinear fe; + FEQuadraticSub fe; PoissonEquation equation (*rhs); QGauss2 quadrature; @@ -481,7 +555,7 @@ void PoissonProblem::run (ParameterHandler &prm) { string o_filename = prm.get ("Output file"); ofstream gnuplot(o_filename.c_str()); fill_data (out); - out.write_gnuplot (gnuplot); + out.write_ucd (gnuplot); gnuplot.close (); // release the lock of the DoF object to @@ -495,4 +569,4 @@ void PoissonProblem::run (ParameterHandler &prm) { -template class PoissonProblem<2>; +template class PoissonProblem<3>; -- 2.39.5