From: Martin Kronbichler Date: Wed, 22 Jun 2016 16:12:22 +0000 (+0200) Subject: Place new points for MappingQ along manifold X-Git-Tag: v8.5.0-rc1~958^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=49f3a01ee6bc39bfd1983d2285caf10171713545;p=dealii.git Place new points for MappingQ along manifold --- diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 6f82615823..0ea518fb13 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -3479,7 +3479,7 @@ namespace { case 2: { - // If two points are passed, these are the two vertices, and + // If two points are passed, these are the two vertices, so // we can only compute degree-1 intermediate points. for (unsigned int i=0; i= 3, ExcImpossibleInDim(spacedim)); + unsigned int m=1; + for ( ; m < n; ++m) + if (m*m*m == n) + break; + // is n a cube number + Assert(m*m*m==n, ExcInternalError()); + + // If m points are passed, two of them are the outer vertices, and we can + // only compute (degree-1)*(degree-1)*(degree-1) intermediate points. + for (unsigned int k=0, c=0; k quadrature(surrounding_points, w); + points[c]=manifold.get_new_point(quadrature); + } + } + } break; + } default: Assert(false, ExcInternalError()); break; @@ -3737,9 +3770,10 @@ add_quad_support_points(const Triangulation<3,3>::cell_iterator &cell, ExcInternalError()); #endif - // if face at boundary, then ask boundary object to return intermediate - // points on it - if (face->at_boundary()) + // if face at boundary or if it belongs to a manifold, then ask the + // boundary/manifold object to return intermediate points on it + if (face->at_boundary() || + dynamic_cast *>(&face->get_manifold()) == NULL) { get_intermediate_points_on_object(face->get_manifold(), line_support_points, face, quad_points); @@ -3877,6 +3911,13 @@ compute_mapping_support_points(const typename Triangulation::cell_ // manifold, otherwise compute them from the points around it if (dim != spacedim) add_quad_support_points(cell, a); + else if (dynamic_cast *>(&cell->get_manifold()) == NULL) + { + std::vector > quad_points (Utilities::fixed_power(polynomial_degree-1)); + get_intermediate_points_on_object(cell->get_manifold(), line_support_points, cell, quad_points); + for (unsigned int i=0; i::cell_ add_quad_support_points (cell, a); // then compute the interior points - add_weighted_interior_points (support_point_weights_on_hex, a); + if (dynamic_cast *>(&cell->get_manifold()) == NULL) + { + std::vector > hex_points (Utilities::fixed_power(polynomial_degree-1)); + get_intermediate_points_on_object(cell->get_manifold(), line_support_points, cell, hex_points); + for (unsigned int i=0; i +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace dealii; + +double f_x(double x_m) +{ + double x = x_m*1000.0; + double y_m = 0.0; + + if (x <= 9.0) + y_m = 0.001*(-28.0 + std::min(28.0, 2.8e1 + 6.775070969851e-3*x*x - 2.124527775800e-3*x*x*x)); + + else if (x > 9.0 && x <= 14.0) + y_m = 0.001*(-28.0 + 2.507355893131e1 + 9.754803562315e-1*x - 1.016116352781e-1*x*x + 1.889794677828e-3*x*x*x); + + else if (x > 14.0 && x <= 20.0) + y_m = 0.001*(-28.0 + 2.579601052357e1 + 8.206693007457e-1*x - 9.055370274339e-2*x*x + 1.626510569859e-3*x*x*x); + + else if (x > 20.0 && x <= 30.0) + y_m = 0.001*(-28.0 + 4.046435022819e1 - 1.379581654948*x + 1.945884504128e-2*x*x - 2.070318932190e-4*x*x*x); + + else if (x > 30.0 && x <= 40.0) + y_m = 0.001*(-28.0 + 1.792461334664e1 + 8.743920332081e-1*x - 5.567361123058e-2*x*x + 6.277731764683e-4*x*x*x); + + else if (x > 40.0 && x <= 54.0) + y_m = 0.001*(-28.0 + std::max(0.0, 5.639011190988e1 - 2.010520359035*x + 1.644919857549e-2*x*x + 2.674976141766e-5*x*x*x)); + + else if (x > 54.0) + y_m = 0.001*(-28.0); + return y_m; +} + +template +class PushForward : public Function +{ +public: + PushForward () : Function(dim, 0.) {} //,component(component) {} + + virtual ~PushForward() {}; + + virtual double value (const Point &p,const unsigned int component = 0) const; + +private: + const double h = 0.028; + + // data from initial block + const double x_max = 4.5*h;//9.0*h; + const double y_max = 2.036*h; + const double z_max = 4.5*h; + + const double y_FoR = h; + + +}; + +template +double PushForward::value(const Point &p,const unsigned int component) const +{ + double result = 0; + + // x component + if (component == 0) + result = p[0]; + + // y component + else if (component == 1) + { + if (p[0] <= x_max/2.0) + result = p[1] + (1 - (p[1] - y_FoR)/y_max)*f_x(p[0]); + + else if (p[0] > x_max/2.0) + result = p[1] + (1 - (p[1] - y_FoR)/y_max)*f_x(x_max - p[0]); + } + + // z component + else if (component == 2) + result = p[2]; + + return result; +} + +template +class PullBack : public Function +{ +public: + PullBack () : Function(dim, 0.) {} //,component(component) {} + + virtual ~PullBack() {}; + + virtual double value (const Point &p,const unsigned int component = 0) const; + +private: + const double h = 0.028; + const double x_max = 4.5*h;//9.0*h; + const double y_max = 2.036*h; + const double z_max = 4.5*h; + + const double y_FoR = h; +}; + +template +double PullBack::value(const Point &p,const unsigned int component) const +{ + double result = 0; + + // x component + if (component == 0) + result = p[0]; + + // y component + else if (component == 1) + { + if (p[0] <= x_max/2.0) + result = (p[1] - f_x(p[0])*(1 + y_FoR/y_max))/(1.0 - f_x(p[0])/y_max); + else if (p[0] > x_max/2.0) + result = (p[1] - f_x(x_max - p[0])*(1 + y_FoR/y_max))/(1.0 - f_x(x_max - p[0])/y_max); + } + + // z component + else if (component == 2) + result = p[2]; + + return result; +} + +template +void create_tria(Triangulation &triangulation, const Manifold &manifold) +{ + const double h = 0.028; + std::vector refinements(dim,1); + refinements[1] = 2; + + Point p1, p2; + p2[0] = 4.5*h;//9.0*h; + p1[1] = h; + p2[1] = 2.018*h;//2.036*h; + if (dim == 3) + { + p1[2] = -2.25*h; + p2[2] = 2.25*h; + } + GridGenerator::hyper_rectangle(triangulation, p1, p2); + + for (typename Triangulation::cell_iterator cell = triangulation.begin(); + cell != triangulation.end(); ++cell) + cell->set_all_manifold_ids(111); + triangulation.set_manifold(111, manifold); + + triangulation.refine_global(1); +} + +template +void test() +{ + deallog << "dim: " << dim << std::endl; + + PushForward push_forward; + PullBack pull_back; + FunctionManifold manifold(push_forward, pull_back); + + Triangulation triangulation; + create_tria(triangulation, manifold); + + for (unsigned mapping_p = 4; mapping_p <= 5; ++mapping_p) + { + deallog << "Mapping degree: " << mapping_p << std::endl; + MappingQ mapping(mapping_p, true); + std::vector > points(Utilities::fixed_power(2)); + for (unsigned int i=0, c=0; i<(dim==2?1:2); ++i) + for (unsigned int j=0; j<2; ++j) + for (unsigned int k=0; k<2; ++k, ++c) + { + points[c][0] = 0.1 + 0.8*k; + points[c][1] = 0.1 + 0.8*j; + if (dim == 3) + points[c][2] = 0.1 + 0.8*i; + } + FE_Nothing dummy; + Quadrature quad(points); + FEValues fe_values(mapping, dummy, quad, update_quadrature_points); + + for (typename Triangulation::active_cell_iterator + cell=triangulation.begin_active(); cell != triangulation.end(); ++cell) + { + fe_values.reinit(cell); + deallog << "Points for cell with first point: " << cell->vertex(0) + << std::endl; + for (unsigned int q=0; q(); + test<3>(); +} diff --git a/tests/fe/mapping_q_manifold_02.output b/tests/fe/mapping_q_manifold_02.output new file mode 100644 index 0000000000..93aba769c5 --- /dev/null +++ b/tests/fe/mapping_q_manifold_02.output @@ -0,0 +1,59 @@ + +DEAL::dim: 2 +DEAL::Mapping degree: 4 +DEAL::Points for cell with first point: 0.0000 0.028000 +DEAL::0.0063000 0.029459 0.056700 0.0018472 0.0063000 0.040854 0.056700 0.018906 +DEAL::Points for cell with first point: 0.063000 0.0000 +DEAL::0.069300 0.0018472 0.11970 0.029459 0.069300 0.018906 0.11970 0.040854 +DEAL::Points for cell with first point: 0.0000 0.042252 +DEAL::0.0063000 0.043703 0.056700 0.023171 0.0063000 0.055097 0.056700 0.040229 +DEAL::Points for cell with first point: 0.063000 0.021252 +DEAL::0.069300 0.023171 0.11970 0.043703 0.069300 0.040229 0.11970 0.055097 +DEAL:: +DEAL::Mapping degree: 5 +DEAL::Points for cell with first point: 0.0000 0.028000 +DEAL::0.0063000 0.029303 0.056700 0.0020929 0.0063000 0.040730 0.056700 0.019101 +DEAL::Points for cell with first point: 0.063000 0.0000 +DEAL::0.069300 0.0020929 0.11970 0.029303 0.069300 0.019101 0.11970 0.040730 +DEAL::Points for cell with first point: 0.0000 0.042252 +DEAL::0.0063000 0.043586 0.056700 0.023353 0.0063000 0.055013 0.056700 0.040361 +DEAL::Points for cell with first point: 0.063000 0.021252 +DEAL::0.069300 0.023353 0.11970 0.043586 0.069300 0.040361 0.11970 0.055013 +DEAL:: +DEAL::dim: 3 +DEAL::Mapping degree: 4 +DEAL::Points for cell with first point: 0.0000 0.028000 -0.063000 +DEAL::0.0063000 0.029459 -0.056700 0.056700 0.0018472 -0.056700 0.0063000 0.040854 -0.056700 0.056700 0.018906 -0.056700 0.0063000 0.029459 -0.0063000 0.056700 0.0018472 -0.0063000 0.0063000 0.040854 -0.0063000 0.056700 0.018906 -0.0063000 +DEAL::Points for cell with first point: 0.063000 0.0000 -0.063000 +DEAL::0.069300 0.0018472 -0.056700 0.11970 0.029459 -0.056700 0.069300 0.018906 -0.056700 0.11970 0.040854 -0.056700 0.069300 0.0018472 -0.0063000 0.11970 0.029459 -0.0063000 0.069300 0.018906 -0.0063000 0.11970 0.040854 -0.0063000 +DEAL::Points for cell with first point: 0.0000 0.042252 -0.063000 +DEAL::0.0063000 0.043703 -0.056700 0.056700 0.023171 -0.056700 0.0063000 0.055097 -0.056700 0.056700 0.040229 -0.056700 0.0063000 0.043703 -0.0063000 0.056700 0.023171 -0.0063000 0.0063000 0.055097 -0.0063000 0.056700 0.040229 -0.0063000 +DEAL::Points for cell with first point: 0.063000 0.021252 -0.063000 +DEAL::0.069300 0.023171 -0.056700 0.11970 0.043703 -0.056700 0.069300 0.040229 -0.056700 0.11970 0.055097 -0.056700 0.069300 0.023171 -0.0063000 0.11970 0.043703 -0.0063000 0.069300 0.040229 -0.0063000 0.11970 0.055097 -0.0063000 +DEAL::Points for cell with first point: 0.0000 0.028000 0.0000 +DEAL::0.0063000 0.029459 0.0063000 0.056700 0.0018472 0.0063000 0.0063000 0.040854 0.0063000 0.056700 0.018906 0.0063000 0.0063000 0.029459 0.056700 0.056700 0.0018472 0.056700 0.0063000 0.040854 0.056700 0.056700 0.018906 0.056700 +DEAL::Points for cell with first point: 0.063000 0.0000 0.0000 +DEAL::0.069300 0.0018472 0.0063000 0.11970 0.029459 0.0063000 0.069300 0.018906 0.0063000 0.11970 0.040854 0.0063000 0.069300 0.0018472 0.056700 0.11970 0.029459 0.056700 0.069300 0.018906 0.056700 0.11970 0.040854 0.056700 +DEAL::Points for cell with first point: 0.0000 0.042252 0.0000 +DEAL::0.0063000 0.043703 0.0063000 0.056700 0.023171 0.0063000 0.0063000 0.055097 0.0063000 0.056700 0.040229 0.0063000 0.0063000 0.043703 0.056700 0.056700 0.023171 0.056700 0.0063000 0.055097 0.056700 0.056700 0.040229 0.056700 +DEAL::Points for cell with first point: 0.063000 0.021252 0.0000 +DEAL::0.069300 0.023171 0.0063000 0.11970 0.043703 0.0063000 0.069300 0.040229 0.0063000 0.11970 0.055097 0.0063000 0.069300 0.023171 0.056700 0.11970 0.043703 0.056700 0.069300 0.040229 0.056700 0.11970 0.055097 0.056700 +DEAL:: +DEAL::Mapping degree: 5 +DEAL::Points for cell with first point: 0.0000 0.028000 -0.063000 +DEAL::0.0063000 0.029303 -0.056700 0.056700 0.0020929 -0.056700 0.0063000 0.040730 -0.056700 0.056700 0.019101 -0.056700 0.0063000 0.029303 -0.0063000 0.056700 0.0020929 -0.0063000 0.0063000 0.040730 -0.0063000 0.056700 0.019101 -0.0063000 +DEAL::Points for cell with first point: 0.063000 0.0000 -0.063000 +DEAL::0.069300 0.0020929 -0.056700 0.11970 0.029303 -0.056700 0.069300 0.019101 -0.056700 0.11970 0.040730 -0.056700 0.069300 0.0020929 -0.0063000 0.11970 0.029303 -0.0063000 0.069300 0.019101 -0.0063000 0.11970 0.040730 -0.0063000 +DEAL::Points for cell with first point: 0.0000 0.042252 -0.063000 +DEAL::0.0063000 0.043586 -0.056700 0.056700 0.023353 -0.056700 0.0063000 0.055013 -0.056700 0.056700 0.040361 -0.056700 0.0063000 0.043586 -0.0063000 0.056700 0.023353 -0.0063000 0.0063000 0.055013 -0.0063000 0.056700 0.040361 -0.0063000 +DEAL::Points for cell with first point: 0.063000 0.021252 -0.063000 +DEAL::0.069300 0.023353 -0.056700 0.11970 0.043586 -0.056700 0.069300 0.040361 -0.056700 0.11970 0.055013 -0.056700 0.069300 0.023353 -0.0063000 0.11970 0.043586 -0.0063000 0.069300 0.040361 -0.0063000 0.11970 0.055013 -0.0063000 +DEAL::Points for cell with first point: 0.0000 0.028000 0.0000 +DEAL::0.0063000 0.029303 0.0063000 0.056700 0.0020929 0.0063000 0.0063000 0.040730 0.0063000 0.056700 0.019101 0.0063000 0.0063000 0.029303 0.056700 0.056700 0.0020929 0.056700 0.0063000 0.040730 0.056700 0.056700 0.019101 0.056700 +DEAL::Points for cell with first point: 0.063000 0.0000 0.0000 +DEAL::0.069300 0.0020929 0.0063000 0.11970 0.029303 0.0063000 0.069300 0.019101 0.0063000 0.11970 0.040730 0.0063000 0.069300 0.0020929 0.056700 0.11970 0.029303 0.056700 0.069300 0.019101 0.056700 0.11970 0.040730 0.056700 +DEAL::Points for cell with first point: 0.0000 0.042252 0.0000 +DEAL::0.0063000 0.043586 0.0063000 0.056700 0.023353 0.0063000 0.0063000 0.055013 0.0063000 0.056700 0.040361 0.0063000 0.0063000 0.043586 0.056700 0.056700 0.023353 0.056700 0.0063000 0.055013 0.056700 0.056700 0.040361 0.056700 +DEAL::Points for cell with first point: 0.063000 0.021252 0.0000 +DEAL::0.069300 0.023353 0.0063000 0.11970 0.043586 0.0063000 0.069300 0.040361 0.0063000 0.11970 0.055013 0.0063000 0.069300 0.023353 0.056700 0.11970 0.043586 0.056700 0.069300 0.040361 0.056700 0.11970 0.055013 0.056700 +DEAL::