From 2c56771dcef98f2a4d7866fe5d36aed835804edd Mon Sep 17 00:00:00 2001 From: carraro Date: Thu, 4 Aug 2011 18:26:57 +0000 Subject: [PATCH] added type 3 git-svn-id: https://svn.dealii.org/trunk@24015 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-47/step-47.cc | 82 +++++++++++++++++++++++------ 1 file changed, 67 insertions(+), 15 deletions(-) diff --git a/deal.II/examples/step-47/step-47.cc b/deal.II/examples/step-47/step-47.cc index b12d04af62..121c0e4dc5 100644 --- a/deal.II/examples/step-47/step-47.cc +++ b/deal.II/examples/step-47/step-47.cc @@ -288,12 +288,14 @@ void LaplaceProblem::assemble_system () endc = dof_handler.end(); std::vector level_set_values; + level_set_values.push_back(-1./2.); level_set_values.push_back(1); - level_set_values.push_back(1); - level_set_values.push_back(1); - level_set_values.push_back(-1); + level_set_values.push_back(-1./2.); + level_set_values.push_back(1.); for (; cell!=endc; ++cell) { + compute_quadrature(quadrature_formula, cell, level_set_values); + assert(0); const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; cell_matrix.reinit (dofs_per_cell, dofs_per_cell); @@ -366,8 +368,23 @@ unsigned int LaplaceProblem::compute_quadrature ( Quadrature plain_qua 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: +--+, ++--, +-+-, -++-, --++, -+-+ + 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; + 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 1; @@ -382,7 +399,6 @@ unsigned int LaplaceProblem::compute_quadrature ( Quadrature plain_qua 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; @@ -391,13 +407,6 @@ unsigned int LaplaceProblem::compute_quadrature ( Quadrature plain_qua 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 @@ -474,14 +483,12 @@ unsigned int LaplaceProblem::compute_quadrature ( Quadrature plain_qua 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; @@ -526,11 +533,56 @@ unsigned int LaplaceProblem::compute_quadrature ( Quadrature plain_qua } - return 2; } - return 100; + // Type three decomposition + // (+--+, ++--, +-+-, -++-, --++, -+-+) + + if (type==3) + { + 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); + + if ( sign_ls[0]==sign_ls[1] && sign_ls[2]==sign_ls[3] ) + { + Pos = 0; + A(1) = level_set_values[0]/((level_set_values[0]-level_set_values[2])); + 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])); + B(0) = level_set_values[2]/((level_set_values[2]-level_set_values[3])); + } + 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; + + return 3; + } + + + + return 0; } -- 2.39.5