From: Joerg Frohne Date: Thu, 26 Sep 2013 15:26:10 +0000 (+0000) Subject: splitting obstacle class into two classes X-Git-Tag: v8.1.0~718 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=84ca4bf8c17510f6ba05cee19a343777abdf669d;p=dealii.git splitting obstacle class into two classes git-svn-id: https://svn.dealii.org/trunk@30959 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc index 604346b30f..23550886df 100644 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@ -494,16 +494,61 @@ namespace Step42 // a function (here a ball). // z_max_domain is the z value of the surface of the work piece template - class Obstacle : public Function + class SphereObstacle : public Function { public: - Obstacle (const std_cxx1x::shared_ptr > &input, - const bool use_read_obstacle, + SphereObstacle (const double z_max_domain) + : + Function(dim), + z_max_domain(z_max_domain) + {} + + virtual + double value (const Point &p, + const unsigned int component = 0) const; + + virtual + void vector_value (const Point &p, + Vector &values) const; + + private: + double z_max_domain; + }; + + + template + double + SphereObstacle::value (const Point &p, + const unsigned int component) const + { + if (component == 0) + return p(0); + 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; + } + + template + void + SphereObstacle::vector_value (const Point &p, + Vector &values) const + { + for (unsigned int c = 0; c < this->n_components; ++c) + values(c) = SphereObstacle::value(p, c); + } + + template + class ChineseObstacle : public Function + { + public: + ChineseObstacle (const std_cxx1x::shared_ptr > &input, const double z_max_domain) : Function(dim), input_obstacle(input), - use_read_obstacle(use_read_obstacle), z_max_domain(z_max_domain) {} @@ -517,14 +562,13 @@ namespace Step42 private: const std_cxx1x::shared_ptr > &input_obstacle; - bool use_read_obstacle; double z_max_domain; }; template double - Obstacle::value (const Point &p, + ChineseObstacle::value (const Point &p, const unsigned int component) const { if (component == 0) @@ -533,29 +577,20 @@ namespace Step42 return p(1); //component==2: - if (use_read_obstacle) - { - 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; - } - else - { - //sphere: - 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; - } + 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; + } template void - Obstacle::vector_value (const Point &p, + ChineseObstacle::vector_value (const Point &p, Vector &values) const { for (unsigned int c = 0; c < this->n_components; ++c) - values(c) = Obstacle::value(p, c); + values(c) = ChineseObstacle::value(p, c); } } @@ -1241,8 +1276,10 @@ namespace Step42 void PlasticityContactProblem::update_solution_and_constraints () { - const EquationData::Obstacle obstacle(input_obstacle, - (obstacle_filename != ""), (base_mesh == "box" ? 1.0 : 0.5)); + const EquationData::SphereObstacle obstacle((base_mesh == "box" ? 1.0 : 0.5)); + +// const EquationData::ChineseObstacle obstacle(input_obstacle, (base_mesh == "box" ? 1.0 : 0.5)); + std::vector vertex_touched(dof_handler.n_dofs(), false); typename DoFHandler::active_cell_iterator cell =