From 59c5738e9f0cf392e0cb30c801da00b273ea46b4 Mon Sep 17 00:00:00 2001 From: bangerth Date: Thu, 4 Aug 2011 18:56:37 +0000 Subject: [PATCH] Move slightly further. git-svn-id: https://svn.dealii.org/trunk@24016 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-47/step-47.cc | 134 +++++++++++++++------------- 1 file changed, 74 insertions(+), 60 deletions(-) diff --git a/deal.II/examples/step-47/step-47.cc b/deal.II/examples/step-47/step-47.cc index 121c0e4dc5..4d43cd3c4b 100644 --- a/deal.II/examples/step-47/step-47.cc +++ b/deal.II/examples/step-47/step-47.cc @@ -64,7 +64,7 @@ class LaplaceProblem private: bool interface_intersects_cell (const typename Triangulation::cell_iterator &cell) const; - unsigned int compute_quadrature(Quadrature plain_quadrature, typename hp::DoFHandler::active_cell_iterator cell, std::vector level_set_values); + unsigned int compute_quadrature(const Quadrature &plain_quadrature, const typename hp::DoFHandler::active_cell_iterator &cell, const std::vector &level_set_values); void append_quadrature(Quadrature plain_quadrature, std::vector > v); void setup_system (); @@ -268,7 +268,7 @@ void LaplaceProblem::assemble_system () { const QGauss quadrature_formula(2); - + FEValues plain_fe_values (fe_collection[0], quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values); @@ -346,15 +346,16 @@ void LaplaceProblem::assemble_system () // the elements on both sides of the discontinuity. The disontinuity line is approximated // by a piece-wise linear interpolation between the intersection of the discontinuity // with the edges of the elements. The vector level_set_values has the values of -// the level set function at the vertices of the elements. From these values can be found -// by linear interpolation the intersections. There are three kind of decomposition that +// the level set function at the vertices of the elements. From these values can be found +// by linear interpolation the intersections. There are three kind of decomposition that // are considered. // Type 1: there is not cut. Type 2: a corner of the element is cut. Type 3: two corners are cut. template -unsigned int LaplaceProblem::compute_quadrature ( Quadrature plain_quadrature, - typename hp::DoFHandler::active_cell_iterator cell, - std::vector level_set_values ) + std::pair > + LaplaceProblem::compute_quadrature (const Quadrature &plain_quadrature, + const typename hp::DoFHandler::active_cell_iterator &cell, + const std::vector &level_set_values ) { unsigned int type = 0; @@ -386,7 +387,8 @@ unsigned int LaplaceProblem::compute_quadrature ( Quadrature plain_qua Point E(0,0); Point F(0,0); - if (type == 1) return 1; + if (type == 1) + return std::pair >(1, plain_quadrature); if (type==2) { @@ -396,7 +398,7 @@ unsigned int LaplaceProblem::compute_quadrature ( Quadrature plain_qua // in type 2 there are 5 subelements Quadrature xfem_quadrature(5*n_q_points); - + std::vector > v(GeometryInfo::vertices_per_cell); if (sign_ls[0]!=sign_ls[1] && sign_ls[0]!=sign_ls[2] && sign_ls[0]!=sign_ls[3]) Pos = 0; @@ -528,12 +530,12 @@ unsigned int LaplaceProblem::compute_quadrature ( Quadrature plain_qua std::cout << std::endl; // create quadrature rule append_quadrature( xfem_quadrature, - vertices ); + vertices ); } } - return 2; + return std::pair >(2, xfem_quadrature); } // Type three decomposition @@ -547,7 +549,7 @@ unsigned int LaplaceProblem::compute_quadrature ( Quadrature plain_qua // in type 2 there are 5 subelements Quadrature xfem_quadrature(5*n_q_points); - + std::vector > v(GeometryInfo::vertices_per_cell); if ( sign_ls[0]==sign_ls[1] && sign_ls[2]==sign_ls[3] ) @@ -565,7 +567,7 @@ unsigned int LaplaceProblem::compute_quadrature ( Quadrature plain_qua else if ( sign_ls[0]==sign_ls[3] && sign_ls[1]==sign_ls[2] ) { std::cout << "Error: the element has two cut lines and this is not allowed" << std::endl; - assert(0); + assert(0); } else { @@ -577,7 +579,9 @@ unsigned int LaplaceProblem::compute_quadrature ( Quadrature plain_qua std::cout << A << std::endl; std::cout << B << std::endl; - return 3; +//TODO: fill xfem_quadrature + + return std::pair >(3, xfem_quadrature); } @@ -587,9 +591,11 @@ unsigned int LaplaceProblem::compute_quadrature ( Quadrature plain_qua } template -void LaplaceProblem::append_quadrature ( Quadrature plain_quadrature, - std::vector > v ) - +void LaplaceProblem::append_quadrature ( const Quadrature &plain_quadrature, + const std::vector > &v, + std::vector > &xfem_points, + std::vector &xfem_weights) + { // Project integration points into sub-elements. // Map F1. @@ -600,69 +606,77 @@ void LaplaceProblem::append_quadrature ( Quadrature plain_quadrature, unsigned int n_v = GeometryInfo::vertices_per_cell; - std::vector > q_points = plain_quadrature.get_points(); + std::vector > q_points = plain_quadrature.get_points(); std::vector > q_transf(q_points.size()); - std::vector W = plain_quadrature.get_weights(); + std::vector W = plain_quadrature.get_weights(); std::vector phi(n_v); - std::vector dphi_dxi(n_v); - std::vector dphi_deta(n_v); + std::vector > grad_phi(n_v); + + const unsigned int n_q_points = plain_quadrature.size(); + + std::vector JxW(n_q_points); - for (unsigned int i=0; i JxW(n_q_points); + default: + Assert (false, ExcNotImplemented()); + } - for ( unsigned int i = 1; i < n_q_points; i++) - { - double dx_dxi = 0.; - double dx_deta = 0.; - double dy_dxi = 0.; - double dy_deta = 0.; - // Calculate Jacobian of transformation - for (unsigned int j = 0; j::vertices_per_cell; j++) - { - dx_dxi += dphi_dxi[j] * v[j](0); - dx_deta += dphi_deta[j] * v[j](0); - dy_dxi += dphi_dxi[j] * v[j](1); - dy_deta += dphi_deta[j] * v[j](1); - } + Tensor<2,dim> jacobian; + + // Calculate Jacobian of transformation + for (unsigned int d=0; d::vertices_per_cell; j++) + jacobian[d][e] += grad_phi[j][d] * v[j](e); + +/* - double detJ = dx_dxi * dy_deta - dx_deta * dy_dxi; - JxW[i] = W[i] * detJ; + for (unsigned int j = 0; j::vertices_per_cell; j++) + { + dx_dxi += dphi_dxi[j] * v[j](0); + dx_deta += dphi_deta[j] * v[j](0); + dy_dxi += dphi_dxi[j] * v[j](1); + dy_deta += dphi_deta[j] * v[j](1); + } +*/ + double detJ = determinant(jacobian); + xfem_weights.push_back (W[i] * detJ); // Map integration points from reference element to subcell of reference elemment - double x = 0.; - double y = 0.; - for (unsigned int j = 0; j::vertices_per_cell; j++) - { - x += v[j](0) * phi[j]; - y += v[j](1) * phi[j]; - } - Point q_prime(x,y); - q_transf.push_back(q_prime); + Point q_prime; + for (unsigned int d=0; d::vertices_per_cell; j++) + q_prime[d] += v[j](d) * phi[j]; + xfem_points.push_back(q_prime); } } -- 2.39.5