From 82aa9549f4e46d48cec4b96c2bbd5bdf6f51cfa1 Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 27 Sep 2013 20:49:39 +0000 Subject: [PATCH] Move around a bit of stuff. More documentation. git-svn-id: https://svn.dealii.org/trunk@30987 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-42/step-42.cc | 370 +++++++++++++++------------- 1 file changed, 194 insertions(+), 176 deletions(-) diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc index 022d0f8891..f8a9b46016 100644 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@ -83,146 +83,6 @@ namespace Step42 { using namespace dealii; -// @sect3{The Input class template} - -// This class has the the only purpose to read in data from a picture file -// stored in pbm ascii format. This data will be bilinearly interpolated and -// provides in this way a function which describes an obstacle. -// -// The data which we read from the file will be stored in a double std::vector -// named obstacle_data. This vector composes the base to calculate a -// piecewise bilinear function as a polynomial interpolation. This will be -// done by obstacle_function (). -// -// In the function run() of the class -// PlasticityContactProblem we create an object of the this class -// which will be used in the class Obstacle to supply the obstacle function in -// update_solution_and_constraints() of the class -// PlasticityContactProblem. -// -// The hx,hy variables denote the spacing between pixels in $x$ -// and $y$ directions. nx,ny are the numbers of pixels in each of -// these directions. get_value() returns the value of the image -// at a given location, interpolated from the adjacent pixel values. - template - class Input - { - public: - Input (const std::string &name); - - double - get_value (const double x, - const double y) const; - - private: - std::vector obstacle_data; - double hx, hy; - int nx, ny; - - double get_pixel_value (const int i, - const int j) const; - }; - - - // The constructor of this class reads in the data that describes - // the obstacle from the given file name. - template - Input::Input(const std::string &name) - : - obstacle_data(0), - hx(0), - hy(0), - nx(0), - ny(0) - { - std::ifstream f(name.c_str()); - - std::string temp; - f >> temp >> nx >> ny; - - AssertThrow(nx > 0 && ny > 0, - ExcMessage ("Invalid file format.")); - - for (int k = 0; k < nx * ny; k++) - { - double val; - f >> val; - obstacle_data.push_back(val); - } - - hx = 1.0 / (nx - 1); - hy = 1.0 / (ny - 1); - - if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) - std::cout << "Read obstacle from file <" << name << ">" << std::endl - << "Resolution of the scanned obstacle picture: " << nx << " x " << ny - << std::endl; - } - - template - double - Input::get_pixel_value (const int i, - const int j) const - { - assert(i >= 0 && i < nx); - assert(j >= 0 && j < ny); - return obstacle_data[nx * (ny - 1 - j) + i]; - } - - - template - double - Input::get_value (const double x, - const double y) const - { - const int ix = std::min(std::max((int) (x / hx), 0), nx-2); - const int iy = std::min(std::max((int) (y / hy), 0), ny-2); - - FullMatrix H(4, 4); - Vector X(4); - Vector b(4); - - double xx, yy; - - xx = ix * hx; - yy = iy * hy; - H(0, 0) = xx; - H(0, 1) = yy; - H(0, 2) = xx * yy; - H(0, 3) = 1.0; - b(0) = get_pixel_value(ix, iy); - - xx = (ix + 1) * hx; - yy = iy * hy; - H(1, 0) = xx; - H(1, 1) = yy; - H(1, 2) = xx * yy; - H(1, 3) = 1.0; - b(1) = get_pixel_value(ix + 1, iy); - - xx = (ix + 1) * hx; - yy = (iy + 1) * hy; - H(2, 0) = xx; - H(2, 1) = yy; - H(2, 2) = xx * yy; - H(2, 3) = 1.0; - b(2) = get_pixel_value(ix + 1, iy + 1); - - xx = ix * hx; - yy = (iy + 1) * hy; - H(3, 0) = xx; - H(3, 1) = yy; - H(3, 2) = xx * yy; - H(3, 3) = 1.0; - b(3) = get_pixel_value(ix, iy + 1); - - H.gauss_jordan(); - H.vmult(X, b); - - return (X(0) * x + X(1) * y + X(2) * x * y + X(3)); - } - - // @sect3{The ConstitutiveLaw class template} // This class provides an interface for a constitutive law, i.e., for the @@ -488,21 +348,20 @@ namespace Step42 values(c) = BoundaryValues::value(p, c); } - // This function is obviously implemented to - // define the obstacle that penetrates our deformable - // body. You can choose between two ways to define - // your obstacle: to read it from a file or to use - // a function (here a ball). - // z_max_domain is the z value of the surface of the work piece + + // @sect4{The SphereObstacle class} + + // The following class is the first of two obstacles that can be + // selected from the input file. It describes a sphere centered + // at position $x=y=0.5, z=z_{\text{surface}}+0.59$ and radius $r=0.6$, + // where $z_{\text{surface}}$ is the vertical position of the (flat) + // surface of the deformable body. The function's value + // returns the location of the obstacle for a given $x,y$ value. template class SphereObstacle : public Function { public: - SphereObstacle (const double z_max_domain) - : - Function(dim), - z_max_domain(z_max_domain) - {} + SphereObstacle (const double z_surface); virtual double value (const Point &p, @@ -513,10 +372,18 @@ namespace Step42 Vector &values) const; private: - double z_max_domain; + const double z_surface; }; + template + SphereObstacle::SphereObstacle (const double z_surface) + : + Function(dim), + z_surface(z_surface) + {} + + template double SphereObstacle::value (const Point &p, @@ -524,14 +391,19 @@ namespace Step42 { if (component == 0) return p(0); - if (component == 1) + else if (component == 1) return p(1); - - //component==2: - return -std::sqrt(0.36 - (p(0) - 0.5) * (p(0) - 0.5) - - (p(1) - 0.5) * (p(1) - 0.5)) + z_max_domain + 0.59; + else if (component == 2) + return (-std::sqrt(0.36 + - (p(0) - 0.5) * (p(0) - 0.5) + - (p(1) - 0.5) * (p(1) - 0.5)) + + z_surface + 0.59); + + Assert (false, ExcNotImplemented()); + return 1e9; // an unreasonable value; ignored in debug mode because of the preceding Assert } + template void SphereObstacle::vector_value (const Point &p, @@ -541,17 +413,151 @@ namespace Step42 values(c) = SphereObstacle::value(p, c); } + // @sect4{The BitmapFile and ChineseObstacle classes} + + // The following two classes describe the obstacle outlined in the introduction, + // i.e., the Chinese character. The first of the two, BitmapFile + // is responsible for reading in data from a picture file + // stored in pbm ascii format. This data will be bilinearly interpolated and + // provides in this way a function which describes an obstacle. + // + // The data which we read from the file will be stored in a double std::vector + // named obstacle_data. This vector composes the base to calculate a + // piecewise bilinear function as a polynomial interpolation. The data we will + // read from a file consists of zeros (white) and ones (black). + // + // The hx,hy variables denote the spacing between pixels in $x$ + // and $y$ directions. nx,ny are the numbers of pixels in each of + // these directions. get_value() returns the value of the image + // at a given location, interpolated from the adjacent pixel values. + template + class BitmapFile + { + public: + BitmapFile(const std::string &name); + + double + get_value(const double x, const double y) const; + + private: + std::vector obstacle_data; + double hx, hy; + int nx, ny; + + double + get_pixel_value(const int i, const int j) const; + }; + + // The constructor of this class reads in the data that describes + // the obstacle from the given file name. + template + BitmapFile::BitmapFile(const std::string &name) + : + obstacle_data(0), + hx(0), + hy(0), + nx(0), + ny(0) + { + std::ifstream f(name.c_str()); + + std::string temp; + f >> temp >> nx >> ny; + + AssertThrow(nx > 0 && ny > 0, ExcMessage("Invalid file format.")); + + for (int k = 0; k < nx * ny; k++) + { + double val; + f >> val; + obstacle_data.push_back(val); + } + + hx = 1.0 / (nx - 1); + hy = 1.0 / (ny - 1); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + std::cout << "Read obstacle from file <" << name << ">" << std::endl + << "Resolution of the scanned obstacle picture: " << nx + << " x " << ny << std::endl; + } + + template + double + BitmapFile::get_pixel_value(const int i, const int j) const + { + assert(i >= 0 && i < nx); + assert(j >= 0 && j < ny); + return obstacle_data[nx * (ny - 1 - j) + i]; + } + + template + double + BitmapFile::get_value(const double x, const double y) const + { + const int ix = std::min(std::max((int) (x / hx), 0), nx - 2); + const int iy = std::min(std::max((int) (y / hy), 0), ny - 2); + + FullMatrix H(4, 4); + Vector X(4); + Vector b(4); + + double xx, yy; + + xx = ix * hx; + yy = iy * hy; + H(0, 0) = xx; + H(0, 1) = yy; + H(0, 2) = xx * yy; + H(0, 3) = 1.0; + b(0) = get_pixel_value(ix, iy); + + xx = (ix + 1) * hx; + yy = iy * hy; + H(1, 0) = xx; + H(1, 1) = yy; + H(1, 2) = xx * yy; + H(1, 3) = 1.0; + b(1) = get_pixel_value(ix + 1, iy); + + xx = (ix + 1) * hx; + yy = (iy + 1) * hy; + H(2, 0) = xx; + H(2, 1) = yy; + H(2, 2) = xx * yy; + H(2, 3) = 1.0; + b(2) = get_pixel_value(ix + 1, iy + 1); + + xx = ix * hx; + yy = (iy + 1) * hy; + H(3, 0) = xx; + H(3, 1) = yy; + H(3, 2) = xx * yy; + H(3, 3) = 1.0; + b(3) = get_pixel_value(ix, iy + 1); + + H.gauss_jordan(); + H.vmult(X, b); + + return (X(0) * x + X(1) * y + X(2) * x * y + X(3)); + } + + // Finally, this is the class that actually uses the class above. It + // has a BitmapFile object as a member that describes the height of the + // obstacle. As mentioned above, the BitmapFile class will provide us + // with a mask, i.e., values that are either zero or one (and, if you + // ask for locations between pixels, values that are interpolated between + // zero and one). This class translates this to heights that are either + // 0.001 below the surface of the deformable body (if the BitmapFile + // class reports a one at this location) or 0.999 above the obstacle (if + // the BitmapFile class reports a zero). The following function should then + // be self-explanatory. template class ChineseObstacle : public Function { public: - ChineseObstacle (const std::string &filename, - const double z_max_domain) - : - Function(dim), - input_obstacle(filename), - z_max_domain(z_max_domain) - {} + ChineseObstacle(const std::string &filename, + const double z_surface); virtual double value (const Point &p, @@ -562,11 +568,21 @@ namespace Step42 Vector &values) const; private: - const Input input_obstacle; - double z_max_domain; + const BitmapFile input_obstacle; + double z_surface; }; + template + ChineseObstacle::ChineseObstacle(const std::string &filename, + const double z_surface) + : + Function(dim), + input_obstacle(filename), + z_surface(z_surface) + {} + + template double ChineseObstacle::value (const Point &p, @@ -576,13 +592,14 @@ namespace Step42 return p(0); if (component == 1) return p(1); - - //component==2: + else if (component==2) + { if (p(0) >= 0.0 && p(0) <= 1.0 && p(1) >= 0.0 && p(1) <= 1.0) - return z_max_domain + 0.999 - input_obstacle.get_value(p(0), p(1)); - else - return 10000.0; + return z_surface + 0.999 - input_obstacle.get_value(p(0), p(1)); + } + Assert (false, ExcNotImplemented()); + return 1e9; // an unreasonable value; ignored in debug mode because of the preceding Assert } template @@ -601,13 +618,13 @@ namespace Step42 // and variables needed to describe // the nonlinear contact problem. It is // close to step-41 but with some additional - // features like: handling hanging nodes, - // a newton method, using Trilinos and p4est + // features like handling hanging nodes, + // a Newton method, using Trilinos and p4est // for parallel distributed computing. // To deal with hanging nodes makes // life a bit more complicated since - // we need an other ConstraintMatrix now. - // We create a newton method for the + // we need another ConstraintMatrix now. + // We create a Newton method for the // active set method for the contact // situation and to handle the nonlinear // operator for the constitutive law. @@ -745,7 +762,9 @@ namespace Step42 throw ExcNotImplemented(); n_cycles = prm.get_integer("number of cycles"); + base_mesh = prm.get("base mesh"); const std::string obstacle_filename = prm.get("obstacle filename"); + if (obstacle_filename != "") obstacle.reset (new EquationData::ChineseObstacle(obstacle_filename, (base_mesh == "box" ? 1.0 : 0.5))); else @@ -757,7 +776,6 @@ namespace Step42 mkdir(output_dir.c_str(), 0777); transfer_solution = prm.get_bool("transfer solution"); - base_mesh = prm.get("base mesh"); pcout << " Using output directory '" << output_dir << "'" << std::endl; pcout << " FE degree " << degree << std::endl; -- 2.39.5