From 4685a72a8edd153ae74e87bc2a8c20ba8da966a0 Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 5 Aug 2011 18:58:05 +0000 Subject: [PATCH] Make it run. Not tested yet. git-svn-id: https://svn.dealii.org/trunk@24022 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-47/step-47.cc | 924 ++++++++++++++++------------ 1 file changed, 526 insertions(+), 398 deletions(-) diff --git a/deal.II/examples/step-47/step-47.cc b/deal.II/examples/step-47/step-47.cc index fadbfeed4c..420b63728a 100644 --- a/deal.II/examples/step-47/step-47.cc +++ b/deal.II/examples/step-47/step-47.cc @@ -53,6 +53,17 @@ using namespace dealii; +double sign (double d) +{ + if (d > 0) + return 1; + else if (d < 0) + return -1; + else + return 0; +} + + template class LaplaceProblem { @@ -63,12 +74,15 @@ class LaplaceProblem void run (); private: + double level_set (const Point &p) const; + Tensor<1,dim> grad_level_set (const Point &p) const; + bool interface_intersects_cell (const typename Triangulation::cell_iterator &cell) const; - std::pair > compute_quadrature(const Quadrature &plain_quadrature, const typename hp::DoFHandler::active_cell_iterator &cell, const std::vector &level_set_values); - void append_quadrature(const Quadrature &plain_quadrature, - const std::vector > &v , - std::vector > &xfem_points, - std::vector &xfem_weights); + std::pair > compute_quadrature(const Quadrature &plain_quadrature, const typename hp::DoFHandler::active_cell_iterator &cell, const std::vector &level_set_values); + void append_quadrature(const Quadrature &plain_quadrature, + const std::vector > &v , + std::vector > &xfem_points, + std::vector &xfem_weights); void setup_system (); void assemble_system (); @@ -167,11 +181,38 @@ LaplaceProblem::~LaplaceProblem () +template +double +LaplaceProblem:: +level_set (const Point &p) const +{ + return p.norm() - 0.5; +} + + + +template +Tensor<1,dim> +LaplaceProblem:: +grad_level_set (const Point &p) const +{ + return p / p.norm(); +} + + + template bool LaplaceProblem:: interface_intersects_cell (const typename Triangulation::cell_iterator &cell) const { + for (unsigned int v=0; v::vertices_per_cell-1; ++v) + if (level_set(cell->vertex(v)) * level_set(cell->vertex(v+1)) < 0) + return true; + + // we get here only if all vertices + // have the same sign, which means + // that the cell is not intersected return false; } @@ -230,9 +271,9 @@ void LaplaceProblem::setup_system () constraints.clear (); - DoFTools::make_hanging_node_constraints (dof_handler, - constraints); - constraints.close (); +//TODO: fix this, it currently crashes + // DoFTools::make_hanging_node_constraints (dof_handler, + // constraints); // now constrain those enriched // DoFs that are on cells that are @@ -252,7 +293,9 @@ void LaplaceProblem::setup_system () // constrain these DoFs for (unsigned int v=0; v::vertices_per_cell; ++v) if (vertex_is_on_intersected_cell[cell->vertex_index(v)] == false) - constraints.add_line (cell->vertex_dof_index(v,1)); + constraints.add_line (cell->vertex_dof_index(v,1,cell->active_fe_index())); + +//TODO: component 1 must satisfy zero boundary conditions constraints.close(); @@ -275,6 +318,9 @@ void LaplaceProblem::assemble_system () FEValues plain_fe_values (fe_collection[0], quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values); + FEValues enriched_fe_values (fe_collection[1], quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); const unsigned int n_q_points = quadrature_formula.size(); @@ -290,30 +336,13 @@ void LaplaceProblem::assemble_system () cell = dof_handler.begin_active(), endc = dof_handler.end(); - std::vector level_set_values; - level_set_values.push_back(1.); - level_set_values.push_back(1.); - level_set_values.push_back(1.); - level_set_values.push_back(-3.); + std::vector level_set_values; + level_set_values.push_back(1.); + level_set_values.push_back(1.); + level_set_values.push_back(1.); + level_set_values.push_back(-3.); for (; cell!=endc; ++cell) { - std::pair > type_and_quadrature = compute_quadrature(quadrature_formula, cell, level_set_values); - - std::cout << "type : " << type_and_quadrature.first << std::endl; - std::vector > points = type_and_quadrature.second.get_points(); - std::vector weights = type_and_quadrature.second.get_weights(); - std::string filename = "points.dat"; - std::ofstream output (filename.c_str()); - output << "#xfem quadrature Points" << std::endl; - for (unsigned int i=0; iget_fe().dofs_per_cell; cell_matrix.reinit (dofs_per_cell, dofs_per_cell); cell_rhs.reinit (dofs_per_cell); @@ -321,24 +350,121 @@ void LaplaceProblem::assemble_system () cell_matrix = 0; cell_rhs = 0; - plain_fe_values.reinit (cell); + if (cell->active_fe_index() == 0) + { + plain_fe_values.reinit (cell); - coefficient.value_list (plain_fe_values.get_quadrature_points(), - coefficient_values); + coefficient_values.resize (plain_fe_values.n_quadrature_points); + coefficient.value_list (plain_fe_values.get_quadrature_points(), + coefficient_values); - for (unsigned int q_point=0; q_pointactive_fe_index() == 1, ExcInternalError()); + + std::auto_ptr > custom_fe_values; + + if (interface_intersects_cell(cell) == false) + enriched_fe_values.reinit (cell); + else + { + std::vector level_set_values (GeometryInfo::vertices_per_cell); + for (unsigned int v=0; v::vertices_per_cell; ++v) + level_set_values[v] = level_set (cell->vertex(v)); + + custom_fe_values + .reset(new FEValues (fe_collection[1], + compute_quadrature(quadrature_formula, cell, + level_set_values).second, + update_values | update_gradients | + update_quadrature_points | update_JxW_values)); + custom_fe_values->reinit (cell); + } + + FEValues &this_fe_values = ((interface_intersects_cell(cell) == false) + ? + enriched_fe_values + : + *custom_fe_values); + + coefficient_values.resize (this_fe_values.n_quadrature_points); + coefficient.value_list (this_fe_values.get_quadrature_points(), + coefficient_values); + + for (unsigned int q_point=0; q_pointget_fe().system_to_component_index(i).first == 0) + { + for (unsigned int j=0; jget_fe().system_to_component_index(j).first == 0) + cell_matrix(i,j) += (coefficient_values[q_point] * + this_fe_values.shape_grad(i,q_point) * + this_fe_values.shape_grad(j,q_point) * + this_fe_values.JxW(q_point)); + else + cell_matrix(i,j) += (coefficient_values[q_point] * + this_fe_values.shape_grad(i,q_point) * + (std::fabs(level_set(this_fe_values.quadrature_point(q_point))) * + this_fe_values.shape_grad(j,q_point) + + + grad_level_set(this_fe_values.quadrature_point(q_point)) * + sign(level_set(this_fe_values.quadrature_point(q_point))) * + this_fe_values.shape_value(j,q_point)) * + this_fe_values.JxW(q_point)); + + cell_rhs(i) += (this_fe_values.shape_value(i,q_point) * + 1.0 * + this_fe_values.JxW(q_point)); + } + else + { + for (unsigned int j=0; jget_fe().system_to_component_index(j).first == 0) + cell_matrix(i,j) += (coefficient_values[q_point] * + (std::fabs(level_set(this_fe_values.quadrature_point(q_point))) * + this_fe_values.shape_grad(i,q_point) + + + grad_level_set(this_fe_values.quadrature_point(q_point)) * + sign(level_set(this_fe_values.quadrature_point(q_point))) * + this_fe_values.shape_value(i,q_point)) * + this_fe_values.shape_grad(j,q_point) * + this_fe_values.JxW(q_point)); + else + cell_matrix(i,j) += (coefficient_values[q_point] * + (std::fabs(level_set(this_fe_values.quadrature_point(q_point))) * + this_fe_values.shape_grad(i,q_point) + + + grad_level_set(this_fe_values.quadrature_point(q_point)) * + sign(level_set(this_fe_values.quadrature_point(q_point))) * + this_fe_values.shape_value(i,q_point)) * + (std::fabs(level_set(this_fe_values.quadrature_point(q_point))) * + this_fe_values.shape_grad(j,q_point) + + + grad_level_set(this_fe_values.quadrature_point(q_point)) * + sign(level_set(this_fe_values.quadrature_point(q_point))) * + this_fe_values.shape_value(j,q_point)) * + this_fe_values.JxW(q_point)); + + cell_rhs(i) += (std::fabs(level_set(this_fe_values.quadrature_point(q_point))) * + this_fe_values.shape_value(i,q_point) * + 1.0 * + this_fe_values.JxW(q_point)); + } + } local_dof_indices.resize (dofs_per_cell); cell->get_dof_indices (local_dof_indices); @@ -368,408 +494,410 @@ void LaplaceProblem::assemble_system () // are considered. // Type 1: there is not cut. Type 2: a corner of the element is cut. Type 3: two corners are cut. - template - std::pair > - LaplaceProblem::compute_quadrature (const Quadrature &plain_quadrature, - const typename hp::DoFHandler::active_cell_iterator &cell, - const std::vector &level_set_values ) +template +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; + unsigned int type = 0; - // find the type of cut - int sign_ls[GeometryInfo::vertices_per_cell]; - for (unsigned int v=0; v::vertices_per_cell; ++v) - { - if (level_set_values[v] > 0) sign_ls[v] = 1; - else if (level_set_values[v] < 0) sign_ls[v] = -1; - else sign_ls[v] = 0; - } + // find the type of cut + int sign_ls[GeometryInfo::vertices_per_cell]; + for (unsigned int v=0; v::vertices_per_cell; ++v) + { + if (level_set_values[v] > 0) sign_ls[v] = 1; + else if (level_set_values[v] < 0) sign_ls[v] = -1; + else sign_ls[v] = 0; + } - // the sign of the level set function at the 4 nodes of the elements can be positive + or negative - - // depending on the sign of the level set function we have the folloing three classes of decomposition - // type 1: ++++, ---- - // type 2: -+++, +-++, ++-+, +++-, +---, -+--, --+-, ---+ - // type 3: +--+, ++--, +-+-, -++-, --++, -+-+ + // the sign of the level set function at the 4 nodes of the elements can be positive + or negative - + // depending on the sign of the level set function we have the folloing three classes of decomposition + // type 1: ++++, ---- + // type 2: -+++, +-++, ++-+, +++-, +---, -+--, --+-, ---+ + // type 3: +--+, ++--, +-+-, -++-, --++, -+-+ - if ( sign_ls[0]==sign_ls[1] & sign_ls[0]==sign_ls[2] & sign_ls[0]==sign_ls[3] ) type =1; - else if ( sign_ls[0]*sign_ls[1]*sign_ls[2]*sign_ls[3] < 0 ) type = 2; - else type = 3; + if ( sign_ls[0]==sign_ls[1] & sign_ls[0]==sign_ls[2] & sign_ls[0]==sign_ls[3] ) type =1; + else if ( sign_ls[0]*sign_ls[1]*sign_ls[2]*sign_ls[3] < 0 ) type = 2; + else type = 3; - unsigned int Pos = 100; + unsigned int Pos = 100; - Point v0(0,0); - Point v1(1,0); - Point v2(0,1); - Point v3(1,1); + Point v0(0,0); + Point v1(1,0); + Point v2(0,1); + Point v3(1,1); - Point A(0,0); - Point B(0,0); - Point C(0,0); - Point D(0,0); - Point E(0,0); - Point F(0,0); + Point A(0,0); + Point B(0,0); + Point C(0,0); + Point D(0,0); + Point E(0,0); + Point F(0,0); - if (type == 1) - return std::pair >(1, plain_quadrature); + if (type == 1) + return std::pair >(1, plain_quadrature); - if (type==2) - { - const unsigned int n_q_points = plain_quadrature.size(); + if (type==2) + { + const unsigned int n_q_points = plain_quadrature.size(); - // loop over all subelements for integration - // in type 2 there are 5 subelements + // loop over all subelements for integration + // in type 2 there are 5 subelements - Quadrature xfem_quadrature(5*n_q_points); + Quadrature xfem_quadrature(5*n_q_points); - std::vector > v(GeometryInfo::vertices_per_cell); + 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; - else if (sign_ls[1]!=sign_ls[0] && sign_ls[1]!=sign_ls[2] && sign_ls[1]!=sign_ls[3]) Pos = 1; - else if (sign_ls[2]!=sign_ls[0] && sign_ls[2]!=sign_ls[1] && sign_ls[2]!=sign_ls[3]) Pos = 2; - else if (sign_ls[3]!=sign_ls[0] && sign_ls[3]!=sign_ls[1] && sign_ls[3]!=sign_ls[2]) Pos = 3; - else assert(0); // error message + if (sign_ls[0]!=sign_ls[1] && sign_ls[0]!=sign_ls[2] && sign_ls[0]!=sign_ls[3]) Pos = 0; + else if (sign_ls[1]!=sign_ls[0] && sign_ls[1]!=sign_ls[2] && sign_ls[1]!=sign_ls[3]) Pos = 1; + else if (sign_ls[2]!=sign_ls[0] && sign_ls[2]!=sign_ls[1] && sign_ls[2]!=sign_ls[3]) Pos = 2; + else if (sign_ls[3]!=sign_ls[0] && sign_ls[3]!=sign_ls[1] && sign_ls[3]!=sign_ls[2]) Pos = 3; + else assert(0); // error message - // Find cut coordinates + // Find cut coordinates - // deal.ii local coordinates + // deal.ii local coordinates - // 2-------3 - // | | - // | | - // | | - // 0-------1 + // 2-------3 + // | | + // | | + // | | + // 0-------1 - if (Pos == 0) - { - A[0] = 1. - level_set_values[1]/(level_set_values[1]-level_set_values[0]); - B[1] = 1. - level_set_values[2]/(level_set_values[2]-level_set_values[0]); - A(1) = 0.; - B(0) = 0.; - C(0) = 0.5*( A(0) + B(0) ); - C(1) = 0.5*( A(1) + B(1) ); - D(0) = 2./3. * C(0); - D(1) = 2./3. * C(1); - E(0) = 0.5*A(0); - E(1) = 0.; - F(0) = 0.; - F(1) = 0.5*B(1); - } - else if (Pos == 1) - { - A[0] = level_set_values[0]/(level_set_values[0]-level_set_values[1]); - B[1] = 1 - level_set_values[3]/(level_set_values[3]-level_set_values[1]); - A(1) = 0.; - B(0) = 1.; - C(0) = 0.5*( A(0) + B(0) ); - C(1) = 0.5*( A(1) + B(1) ); - D(0) = 1./3. + 2./3. * C(0); - D(1) = 2./3. * C(1); - E(0) = 0.5*(1 + A(0)); - E(1) = 0.; - F(0) = 1.; - F(1) = 0.5*B(1); - } - else if (Pos == 2) - { - A[0] = 1 - level_set_values[3]/(level_set_values[3]-level_set_values[2]); - B[1] = level_set_values[0]/(level_set_values[0]-level_set_values[2]); - A(1) = 1.; - B(0) = 0.; - C(0) = 0.5*( A(0) + B(0) ); - C(1) = 0.5*( A(1) + B(1) ); - D(0) = 2./3. * C(0); - D(1) = 1./3. + 2./3. * C(1); - E(0) = 0.5* A(0); - E(1) = 1.; - F(0) = 0.; - F(1) = 0.5*( 1. + B(1) ); - } - else if (Pos == 3) - { - A[0] = level_set_values[2]/(level_set_values[2]-level_set_values[3]); - B[1] = level_set_values[1]/(level_set_values[1]-level_set_values[3]); - A(1) = 1.; - B(0) = 1.; - C(0) = 0.5*( A(0) + B(0) ); - C(1) = 0.5*( A(1) + B(1) ); - D(0) = 1./3. + 2./3. * C(0); - D(1) = 1./3. + 2./3. * C(1); - E(0) = 0.5*( 1. + A(0) ); - E(1) = 1.; - F(0) = 1.; - F(1) = 0.5*( 1. + B(1) ); - } - - //std::cout << A << std::endl; - //std::cout << B << std::endl; - //std::cout << C << std::endl; - //std::cout << D << std::endl; - //std::cout << E << std::endl; - //std::cout << F << std::endl; - - std::string filename = "vertices.dat"; - std::ofstream output (filename.c_str()); - output << "#vertices of xfem subcells" << std::endl; - output << v0(0) << " " << v0(1) << std::endl; - output << v1(0) << " " << v1(1) << std::endl; - output << v3(0) << " " << v3(1) << std::endl; - output << v2(0) << " " << v2(1) << std::endl; - output << std::endl; - output << A(0) << " " << A(1) << std::endl; - output << B(0) << " " << B(1) << std::endl; - output << std::endl; - output << C(0) << " " << C(1) << std::endl; - output << D(0) << " " << D(1) << std::endl; - output << std::endl; - output << D(0) << " " << D(1) << std::endl; - output << E(0) << " " << E(1) << std::endl; - output << std::endl; - output << D(0) << " " << D(1) << std::endl; - output << F(0) << " " << F(1) << std::endl; - output << std::endl; - - if (Pos==0) - output << v3(0) << " " << v3(1) << std::endl; - else if (Pos==1) - output << v2(0) << " " << v2(1) << std::endl; - else if (Pos==2) - output << v1(0) << " " << v1(1) << std::endl; - else if (Pos==3) - output << v0(0) << " " << v0(1) << std::endl; - output << C(0) << " " << C(1) << std::endl; - - Point subcell_vertices[10]; - subcell_vertices[0] = v0; - subcell_vertices[1] = v1; - subcell_vertices[2] = v2; - subcell_vertices[3] = v3; - subcell_vertices[4] = A; - subcell_vertices[5] = B; - subcell_vertices[6] = C; - subcell_vertices[7] = D; - subcell_vertices[8] = E; - subcell_vertices[9] = F; - - std::vector > xfem_points; - std::vector xfem_weights; - - // lookup table for the decomposition - - if (dim==2) + if (Pos == 0) + { + A[0] = 1. - level_set_values[1]/(level_set_values[1]-level_set_values[0]); + B[1] = 1. - level_set_values[2]/(level_set_values[2]-level_set_values[0]); + A(1) = 0.; + B(0) = 0.; + C(0) = 0.5*( A(0) + B(0) ); + C(1) = 0.5*( A(1) + B(1) ); + D(0) = 2./3. * C(0); + D(1) = 2./3. * C(1); + E(0) = 0.5*A(0); + E(1) = 0.; + F(0) = 0.; + F(1) = 0.5*B(1); + } + else if (Pos == 1) + { + A[0] = level_set_values[0]/(level_set_values[0]-level_set_values[1]); + B[1] = 1 - level_set_values[3]/(level_set_values[3]-level_set_values[1]); + A(1) = 0.; + B(0) = 1.; + C(0) = 0.5*( A(0) + B(0) ); + C(1) = 0.5*( A(1) + B(1) ); + D(0) = 1./3. + 2./3. * C(0); + D(1) = 2./3. * C(1); + E(0) = 0.5*(1 + A(0)); + E(1) = 0.; + F(0) = 1.; + F(1) = 0.5*B(1); + } + else if (Pos == 2) + { + A[0] = 1 - level_set_values[3]/(level_set_values[3]-level_set_values[2]); + B[1] = level_set_values[0]/(level_set_values[0]-level_set_values[2]); + A(1) = 1.; + B(0) = 0.; + C(0) = 0.5*( A(0) + B(0) ); + C(1) = 0.5*( A(1) + B(1) ); + D(0) = 2./3. * C(0); + D(1) = 1./3. + 2./3. * C(1); + E(0) = 0.5* A(0); + E(1) = 1.; + F(0) = 0.; + F(1) = 0.5*( 1. + B(1) ); + } + else if (Pos == 3) + { + A[0] = level_set_values[2]/(level_set_values[2]-level_set_values[3]); + B[1] = level_set_values[1]/(level_set_values[1]-level_set_values[3]); + A(1) = 1.; + B(0) = 1.; + C(0) = 0.5*( A(0) + B(0) ); + C(1) = 0.5*( A(1) + B(1) ); + D(0) = 1./3. + 2./3. * C(0); + D(1) = 1./3. + 2./3. * C(1); + E(0) = 0.5*( 1. + A(0) ); + E(1) = 1.; + F(0) = 1.; + F(1) = 0.5*( 1. + B(1) ); + } + + //std::cout << A << std::endl; + //std::cout << B << std::endl; + //std::cout << C << std::endl; + //std::cout << D << std::endl; + //std::cout << E << std::endl; + //std::cout << F << std::endl; + + std::string filename = "vertices.dat"; + std::ofstream output (filename.c_str()); + output << "#vertices of xfem subcells" << std::endl; + output << v0(0) << " " << v0(1) << std::endl; + output << v1(0) << " " << v1(1) << std::endl; + output << v3(0) << " " << v3(1) << std::endl; + output << v2(0) << " " << v2(1) << std::endl; + output << std::endl; + output << A(0) << " " << A(1) << std::endl; + output << B(0) << " " << B(1) << std::endl; + output << std::endl; + output << C(0) << " " << C(1) << std::endl; + output << D(0) << " " << D(1) << std::endl; + output << std::endl; + output << D(0) << " " << D(1) << std::endl; + output << E(0) << " " << E(1) << std::endl; + output << std::endl; + output << D(0) << " " << D(1) << std::endl; + output << F(0) << " " << F(1) << std::endl; + output << std::endl; + + if (Pos==0) + output << v3(0) << " " << v3(1) << std::endl; + else if (Pos==1) + output << v2(0) << " " << v2(1) << std::endl; + else if (Pos==2) + output << v1(0) << " " << v1(1) << std::endl; + else if (Pos==3) + output << v0(0) << " " << v0(1) << std::endl; + output << C(0) << " " << C(1) << std::endl; + + Point subcell_vertices[10]; + subcell_vertices[0] = v0; + subcell_vertices[1] = v1; + subcell_vertices[2] = v2; + subcell_vertices[3] = v3; + subcell_vertices[4] = A; + subcell_vertices[5] = B; + subcell_vertices[6] = C; + subcell_vertices[7] = D; + subcell_vertices[8] = E; + subcell_vertices[9] = F; + + std::vector > xfem_points; + std::vector xfem_weights; + + // lookup table for the decomposition + + if (dim==2) + { + unsigned int subcell_v_indices[4][5][4] = { + {{0,8,9,7}, {9,7,5,6}, {8,4,7,6}, {5,6,2,3}, {6,4,3,1}}, + {{8,1,7,9}, {4,8,6,7}, {6,7,5,9}, {0,4,2,6}, {2,6,3,5}}, + {{9,7,2,8}, {5,6,9,7}, {6,4,7,8}, {0,1,5,6}, {6,1,4,3}}, + {{7,9,8,3}, {4,6,8,7}, {6,5,7,9}, {0,6,2,4}, {0,1,6,5}} + }; + + std::cout << "Pos : " << Pos << std::endl; + for (unsigned int subcell = 0; subcell<5; subcell++) + { + //std::cout << "subcell : " << subcell << std::endl; + std::vector > vertices; + for (unsigned int i=0; i<4; i++) { - unsigned int subcell_v_indices[4][5][4] = { - {{0,8,9,7}, {9,7,5,6}, {8,4,7,6}, {5,6,2,3}, {6,4,3,1}}, - {{8,1,7,9}, {4,8,6,7}, {6,7,5,9}, {0,4,2,6}, {2,6,3,5}}, - {{9,7,2,8}, {5,6,9,7}, {6,4,7,8}, {0,1,5,6}, {6,1,4,3}}, - {{7,9,8,3}, {4,6,8,7}, {6,5,7,9}, {0,6,2,4}, {0,1,6,5}} - }; - - std::cout << "Pos : " << Pos << std::endl; - for (unsigned int subcell = 0; subcell<5; subcell++) - { - //std::cout << "subcell : " << subcell << std::endl; - std::vector > vertices; - for (unsigned int i=0; i<4; i++) - { - vertices.push_back( subcell_vertices[subcell_v_indices[Pos][subcell][i]] ); - //std::cout << "i : " << i << std::endl; - //std::cout << "subcell v : " << subcell_v_indices[Pos][subcell][i] << std::endl; - //std::cout << vertices[i](0) << " " << vertices[i](1) << std::endl; - } - //std::cout << std::endl; - // create quadrature rule - append_quadrature( plain_quadrature, - vertices, - xfem_points, - xfem_weights); - //initialize xfem_quadrature with quadrature points of all subelements - xfem_quadrature.initialize(xfem_points, xfem_weights); - } + vertices.push_back( subcell_vertices[subcell_v_indices[Pos][subcell][i]] ); + //std::cout << "i : " << i << std::endl; + //std::cout << "subcell v : " << subcell_v_indices[Pos][subcell][i] << std::endl; + //std::cout << vertices[i](0) << " " << vertices[i](1) << std::endl; } - - return std::pair >(2, xfem_quadrature); + //std::cout << std::endl; + // create quadrature rule + append_quadrature( plain_quadrature, + vertices, + xfem_points, + xfem_weights); + //initialize xfem_quadrature with quadrature points of all subelements + xfem_quadrature.initialize(xfem_points, xfem_weights); + } } - // Type three decomposition - // (+--+, ++--, +-+-, -++-, --++, -+-+) + Assert (xfem_quadrature.size() == 20, ExcInternalError()); + return std::pair >(2, xfem_quadrature); + } - if (type==3) - { - const unsigned int n_q_points = plain_quadrature.size(); + // Type three decomposition + // (+--+, ++--, +-+-, -++-, --++, -+-+) - // loop over all subelements for integration - // in type 2 there are 5 subelements + if (type==3) + { + const unsigned int n_q_points = plain_quadrature.size(); - Quadrature xfem_quadrature(5*n_q_points); + // loop over all subelements for integration + // in type 2 there are 5 subelements - std::vector > v(GeometryInfo::vertices_per_cell); + Quadrature xfem_quadrature(5*n_q_points); - if ( sign_ls[0]==sign_ls[1] && sign_ls[2]==sign_ls[3] ) - { - Pos = 0; - A(0) = 0.; - A(1) = level_set_values[0]/((level_set_values[0]-level_set_values[2])); - B(0) = 1.; - B(1) = level_set_values[1]/((level_set_values[1]-level_set_values[3])); - } - else if ( sign_ls[0]==sign_ls[2] && sign_ls[1]==sign_ls[3] ) - { - Pos = 1; - A(0) = level_set_values[0]/((level_set_values[0]-level_set_values[1])); - A(1) = 0.; - B(0) = level_set_values[2]/((level_set_values[2]-level_set_values[3])); - B(1) = 1.; - } - 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); - } - else - { - std::cout << "Error: the level set function has not the right values" << std::endl; - assert(0); - } + std::vector > v(GeometryInfo::vertices_per_cell); - //std::cout << "Pos " << Pos << std::endl; - //std::cout << A << std::endl; - //std::cout << B << std::endl; - std::string filename = "vertices.dat"; - std::ofstream output (filename.c_str()); - output << "#vertices of xfem subcells" << std::endl; - output << A(0) << " " << A(1) << std::endl; - output << B(0) << " " << B(1) << std::endl; - - //fill xfem_quadrature - Point subcell_vertices[6]; - subcell_vertices[0] = v0; - subcell_vertices[1] = v1; - subcell_vertices[2] = v2; - subcell_vertices[3] = v3; - subcell_vertices[4] = A; - subcell_vertices[5] = B; - - std::vector > xfem_points; - std::vector xfem_weights; - - if (dim==2) + if ( sign_ls[0]==sign_ls[1] && sign_ls[2]==sign_ls[3] ) + { + Pos = 0; + A(0) = 0.; + A(1) = level_set_values[0]/((level_set_values[0]-level_set_values[2])); + B(0) = 1.; + B(1) = level_set_values[1]/((level_set_values[1]-level_set_values[3])); + } + else if ( sign_ls[0]==sign_ls[2] && sign_ls[1]==sign_ls[3] ) + { + Pos = 1; + A(0) = level_set_values[0]/((level_set_values[0]-level_set_values[1])); + A(1) = 0.; + B(0) = level_set_values[2]/((level_set_values[2]-level_set_values[3])); + B(1) = 1.; + } + 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); + } + else + { + std::cout << "Error: the level set function has not the right values" << std::endl; + assert(0); + } + + //std::cout << "Pos " << Pos << std::endl; + //std::cout << A << std::endl; + //std::cout << B << std::endl; + std::string filename = "vertices.dat"; + std::ofstream output (filename.c_str()); + output << "#vertices of xfem subcells" << std::endl; + output << A(0) << " " << A(1) << std::endl; + output << B(0) << " " << B(1) << std::endl; + + //fill xfem_quadrature + Point subcell_vertices[6]; + subcell_vertices[0] = v0; + subcell_vertices[1] = v1; + subcell_vertices[2] = v2; + subcell_vertices[3] = v3; + subcell_vertices[4] = A; + subcell_vertices[5] = B; + + std::vector > xfem_points; + std::vector xfem_weights; + + if (dim==2) + { + unsigned int subcell_v_indices[2][2][4] = { + {{0,1,4,5}, {4,5,2,3}}, + {{0,4,2,5}, {4,1,5,3}} + }; + + //std::cout << "Pos : " << Pos << std::endl; + for (unsigned int subcell = 0; subcell<2; subcell++) + { + //std::cout << "subcell : " << subcell << std::endl; + std::vector > vertices; + for (unsigned int i=0; i<4; i++) { - unsigned int subcell_v_indices[2][2][4] = { - {{0,1,4,5}, {4,5,2,3}}, - {{0,4,2,5}, {4,1,5,3}} - }; - - //std::cout << "Pos : " << Pos << std::endl; - for (unsigned int subcell = 0; subcell<2; subcell++) - { - //std::cout << "subcell : " << subcell << std::endl; - std::vector > vertices; - for (unsigned int i=0; i<4; i++) - { - vertices.push_back( subcell_vertices[subcell_v_indices[Pos][subcell][i]] ); - //std::cout << "i : " << i << std::endl; - //std::cout << "subcell v : " << subcell_v_indices[Pos][subcell][i] << std::endl; - //std::cout << vertices[i](0) << " " << vertices[i](1) << std::endl; - } - //std::cout << std::endl; - // create quadrature rule - append_quadrature( plain_quadrature, - vertices, - xfem_points, - xfem_weights); - //initialize xfem_quadrature with quadrature points of all subelements - xfem_quadrature.initialize(xfem_points, xfem_weights); - } + vertices.push_back( subcell_vertices[subcell_v_indices[Pos][subcell][i]] ); + //std::cout << "i : " << i << std::endl; + //std::cout << "subcell v : " << subcell_v_indices[Pos][subcell][i] << std::endl; + //std::cout << vertices[i](0) << " " << vertices[i](1) << std::endl; } - return std::pair >(3, xfem_quadrature); + //std::cout << std::endl; + // create quadrature rule + append_quadrature( plain_quadrature, + vertices, + xfem_points, + xfem_weights); + //initialize xfem_quadrature with quadrature points of all subelements + xfem_quadrature.initialize(xfem_points, xfem_weights); + } } + Assert (xfem_quadrature.size() == 8, ExcInternalError()); + return std::pair >(3, xfem_quadrature); + } - return std::pair >(0, plain_quadrature);; + return std::pair >(0, plain_quadrature);; } - template +template 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. - // This maps quadrature points from a reference element to a subelement of a reference element. - // To implement the action of this map the coordinates of the subelements have been calculated (A(0)...F(0),A(1)...F(1)) - // the coordinates of the quadrature points are given by the bi-linear map defined by the form functions - // $x^\prime_i = \sum_j v^\prime \phi_j(x^hat_i)$, where the $\phi_j$ are the shape functions of the FEQ. + // Project integration points into sub-elements. + // This maps quadrature points from a reference element to a subelement of a reference element. + // To implement the action of this map the coordinates of the subelements have been calculated (A(0)...F(0),A(1)...F(1)) + // the coordinates of the quadrature points are given by the bi-linear map defined by the form functions + // $x^\prime_i = \sum_j v^\prime \phi_j(x^hat_i)$, where the $\phi_j$ are the shape functions of the FEQ. - unsigned int n_v = GeometryInfo::vertices_per_cell; + unsigned int n_v = GeometryInfo::vertices_per_cell; - std::vector > q_points = plain_quadrature.get_points(); - std::vector > q_transf(q_points.size()); - std::vector W = plain_quadrature.get_weights(); - std::vector phi(n_v); - std::vector > grad_phi(n_v); + std::vector > q_points = plain_quadrature.get_points(); + std::vector > q_transf(q_points.size()); + std::vector W = plain_quadrature.get_weights(); + std::vector phi(n_v); + std::vector > grad_phi(n_v); - const unsigned int n_q_points = plain_quadrature.size(); + const unsigned int n_q_points = plain_quadrature.size(); - std::vector JxW(n_q_points); + std::vector JxW(n_q_points); - for ( unsigned int i = 0; i < n_q_points; i++) + for ( unsigned int i = 0; i < n_q_points; i++) + { + switch (dim) { - switch (dim) - { - case 2: - { - double xi = q_points[i](0); - double eta = q_points[i](1); - - // Define shape functions on reference element - // we consider a bi-linear mapping - phi[0] = (1. - xi) * (1. - eta); - phi[1] = xi * (1. - eta); - phi[2] = (1. - xi) * eta; - phi[3] = xi * eta; - - grad_phi[0][0] = (-1. + eta); - grad_phi[1][0] = (1. - eta); - grad_phi[2][0] = -eta; - grad_phi[3][0] = eta; - - grad_phi[0][1] = (-1. + xi); - grad_phi[1][1] = -xi; - grad_phi[2][1] = 1-xi; - grad_phi[3][1] = xi; - - break; - } + case 2: + { + double xi = q_points[i](0); + double eta = q_points[i](1); + + // Define shape functions on reference element + // we consider a bi-linear mapping + phi[0] = (1. - xi) * (1. - eta); + phi[1] = xi * (1. - eta); + phi[2] = (1. - xi) * eta; + phi[3] = xi * eta; + + grad_phi[0][0] = (-1. + eta); + grad_phi[1][0] = (1. - eta); + grad_phi[2][0] = -eta; + grad_phi[3][0] = eta; + + grad_phi[0][1] = (-1. + xi); + grad_phi[1][1] = -xi; + grad_phi[2][1] = 1-xi; + grad_phi[3][1] = xi; - default: - Assert (false, ExcNotImplemented()); - } + break; + } + + default: + Assert (false, ExcNotImplemented()); + } - Tensor<2,dim> jacobian; + Tensor<2,dim> jacobian; - // Calculate Jacobian of transformation - for (unsigned int d=0; d::vertices_per_cell; j++) - { + // Calculate Jacobian of transformation + for (unsigned int d=0; d::vertices_per_cell; j++) + { jacobian[d][e] += grad_phi[j][e] * v[j](d); - } - } - - double detJ = determinant(jacobian); - xfem_weights.push_back (W[i] * detJ); - - // Map integration points from reference element to subcell of reference element - 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); - } + } + } + + double detJ = determinant(jacobian); + xfem_weights.push_back (W[i] * detJ); + + // Map integration points from reference element to subcell of reference element + 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); + } } @@ -848,7 +976,7 @@ void LaplaceProblem::run () static const HyperBallBoundary boundary; triangulation.set_boundary (0, boundary); - triangulation.refine_global (1); + triangulation.refine_global (3); } else refine_grid (); -- 2.39.5