From: carraro Date: Thu, 4 Aug 2011 15:32:09 +0000 (+0000) Subject: Step 47 XFEM X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8de2f89a757275fc314ecffe2ada648fac5da994;p=dealii-svn.git Step 47 XFEM append_quadrature implemented git-svn-id: https://svn.dealii.org/trunk@24014 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-3/step-3.cc b/deal.II/examples/step-3/step-3.cc index b3bb984701..10612cf5eb 100644 --- a/deal.II/examples/step-3/step-3.cc +++ b/deal.II/examples/step-3/step-3.cc @@ -210,7 +210,7 @@ void Step3::make_grid () // has 32 times 32 cells, for a // total of 1024. GridGenerator::hyper_cube (triangulation, -1, 1); - triangulation.refine_global (5); + triangulation.refine_global (0); // Unsure that 1024 is the correct number? // Let's see: n_active_cells return the // number of active cells: diff --git a/deal.II/examples/step-47/Makefile b/deal.II/examples/step-47/Makefile index 575312d3ad..926ea1ca23 100644 --- a/deal.II/examples/step-47/Makefile +++ b/deal.II/examples/step-47/Makefile @@ -105,7 +105,7 @@ clean: # and write the string indicating what we do instead. ./%.g.$(OBJEXT) : @echo "==============debug========= $( $@" - @$(CXX) $(CXXFLAGS.g) -c $< -o $@ + @$(CXX) -w $(CXXFLAGS.g) -c $< -o $@ ./%.$(OBJEXT) : @echo "==============optimized===== $( $@" @$(CXX) $(CXXFLAGS.o) -c $< -o $@ diff --git a/deal.II/examples/step-47/step-47.cc b/deal.II/examples/step-47/step-47.cc index 3ca5026288..b12d04af62 100644 --- a/deal.II/examples/step-47/step-47.cc +++ b/deal.II/examples/step-47/step-47.cc @@ -64,6 +64,8 @@ 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); + void append_quadrature(Quadrature plain_quadrature, std::vector > v); void setup_system (); void assemble_system (); @@ -261,12 +263,12 @@ void LaplaceProblem::setup_system () system_matrix.reinit (sparsity_pattern); } - template void LaplaceProblem::assemble_system () { - const QGauss quadrature_formula(3); + const QGauss quadrature_formula(2); + FEValues plain_fe_values (fe_collection[0], quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values); @@ -284,8 +286,15 @@ void LaplaceProblem::assemble_system () typename hp::DoFHandler::active_cell_iterator 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(-1); for (; cell!=endc; ++cell) { + const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; cell_matrix.reinit (dofs_per_cell, dofs_per_cell); cell_rhs.reinit (dofs_per_cell); @@ -330,8 +339,281 @@ void LaplaceProblem::assemble_system () system_rhs); } +// To integrate the enriched elements we have to find the geometrical decomposition +// of the original element in subelements. The subelements are used to integrate +// 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 +// 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 ) +{ + + 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; + } + + 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; + + if (type == 1) return 1; + + 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 + + Quadrature xfem_quadrature(5*n_q_points); + + std::vector > v(GeometryInfo::vertices_per_cell); + + unsigned int Pos = 100; + 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 + + std::cout << "Pos " << Pos << std::endl; + + Point A(0,0); + Point B(0,0); + Point C(0,0); + Point D(0,0); + Point E(0,0); + Point F(0,0); + + // Find cut coordinates + + // deal.ii local coordinates + + // 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) = 2./3. * C(0); + D(1) = 2./3. * C(1); + E(0) = 0.5*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) = 0.; + 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) ); + } + + Point v0(0,0); + Point v1(1,0); + Point v2(0,1); + Point v3(1,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; + */ + + 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; + + // 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}, {4,1,6,3}}, + {{8,1,7,9}, {4,8,6,7}, {7,9,6,8}, {0,4,2,6}, {2,6,3,5}}, + {{7,9,8,3}, {4,6,8,7}, {6,5,7,9}, {0,6,2,4}, {0,1,6,5}}, + {{9,7,2,8}, {5,6,9,7}, {6,4,7,8}, {0,1,5,6}, {6,1,4,3}} + }; + + for (unsigned int subcell = 0; subcell<5; subcell++) + { + std::vector > vertices; + for (unsigned int i=0; i<4; i++) + { + vertices.push_back( subcell_vertices[subcell_v_indices[Pos][subcell][i]] ); + //std::cout << "Pos : " << Pos << std::endl; + //std::cout << "subcell : " << subcell << std::endl; + //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( xfem_quadrature, + vertices ); + } + + } + + + return 2; + } + + return 100; + +} + + template +void LaplaceProblem::append_quadrature ( Quadrature plain_quadrature, + std::vector > v ) + +{ + // Project integration points into sub-elements. + // Map F1. + // The map F1 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; + + 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 dphi_dxi(n_v); + std::vector dphi_deta(n_v); + + for (unsigned int i=0; i JxW(n_q_points); + + 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); + } + + double detJ = dx_dxi * dy_deta - dx_deta * dy_dxi; + JxW[i] = 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); + } + +} template