From: Wolfgang Bangerth Date: Wed, 23 May 2018 22:35:23 +0000 (+0800) Subject: Remove the never finished step-47. X-Git-Tag: v9.1.0-rc1~1106^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c84f5cf1b638b52f518e081617e9084c325d919f;p=dealii.git Remove the never finished step-47. We started this program in 2012, but it was never finished and there is no plan to ever do so. A paper on the topic using deal.II has also been published in the meantime: Carraro, Wetterauer: On the implementation of the eXtended Finite Element Method (XFEM) for interface problems Archive of Numerical Software, vol. 4, no. 2. --- diff --git a/examples/step-47/CMakeLists.txt b/examples/step-47/CMakeLists.txt deleted file mode 100644 index 4b1d0dc89f..0000000000 --- a/examples/step-47/CMakeLists.txt +++ /dev/null @@ -1,39 +0,0 @@ -## -# CMake script for the step-47 tutorial program: -## - -# Set the name of the project and target: -SET(TARGET "step-47") - -# Declare all source files the target consists of. Here, this is only -# the one step-X.cc file, but as you expand your project you may wish -# to add other source files as well. If your project becomes much larger, -# you may want to either replace the following statement by something like -# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc") -# FILE(GLOB_RECURSE TARGET_INC "include/*.h") -# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC}) -# or switch altogether to the large project CMakeLists.txt file discussed -# in the "CMake in user projects" page accessible from the "User info" -# page of the documentation. -SET(TARGET_SRC - ${TARGET}.cc - ) - -# Usually, you will not need to modify anything beyond this point... - -CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12) - -FIND_PACKAGE(deal.II 9.1.0 QUIET - HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} - ) -IF(NOT ${deal.II_FOUND}) - MESSAGE(FATAL_ERROR "\n" - "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" - "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" - "or set an environment variable \"DEAL_II_DIR\" that contains this path." - ) -ENDIF() - -DEAL_II_INITIALIZE_CACHED_VARIABLES() -PROJECT(${TARGET}) -DEAL_II_INVOKE_AUTOPILOT() diff --git a/examples/step-47/doc/builds-on b/examples/step-47/doc/builds-on deleted file mode 100644 index 9fdd726254..0000000000 --- a/examples/step-47/doc/builds-on +++ /dev/null @@ -1 +0,0 @@ -step-5 diff --git a/examples/step-47/doc/intro.dox b/examples/step-47/doc/intro.dox deleted file mode 100644 index 1518606e43..0000000000 --- a/examples/step-47/doc/intro.dox +++ /dev/null @@ -1,3 +0,0 @@ - -

Introduction

- diff --git a/examples/step-47/doc/kind b/examples/step-47/doc/kind deleted file mode 100644 index 6816e9090f..0000000000 --- a/examples/step-47/doc/kind +++ /dev/null @@ -1 +0,0 @@ -unfinished diff --git a/examples/step-47/doc/points.dat b/examples/step-47/doc/points.dat deleted file mode 100644 index 01ce16ad3c..0000000000 --- a/examples/step-47/doc/points.dat +++ /dev/null @@ -1,21 +0,0 @@ -#xfem quadrature Points -0.781998 0.781998 -0.941587 0.72508 -0.72508 0.941587 -0.926335 0.926335 -0.397329 0.926335 -0.583333 0.72508 -0.583333 0.941587 -0.686004 0.781998 -0.72508 0.583333 -0.926335 0.397329 -0.781998 0.686004 -0.941587 0.583333 -0.115331 0.315492 -0.430422 0.60008 -0.069578 0.816587 -0.259669 0.892842 -0.315492 0.115331 -0.816587 0.069578 -0.60008 0.430422 -0.892842 0.259669 diff --git a/examples/step-47/doc/results.dox b/examples/step-47/doc/results.dox deleted file mode 100644 index f4c6feefb5..0000000000 --- a/examples/step-47/doc/results.dox +++ /dev/null @@ -1 +0,0 @@ -

Results

diff --git a/examples/step-47/doc/tooltip b/examples/step-47/doc/tooltip deleted file mode 100644 index f2003d7a64..0000000000 --- a/examples/step-47/doc/tooltip +++ /dev/null @@ -1 +0,0 @@ -X-FEM diff --git a/examples/step-47/doc/vertices.dat b/examples/step-47/doc/vertices.dat deleted file mode 100644 index 8d1ef87323..0000000000 --- a/examples/step-47/doc/vertices.dat +++ /dev/null @@ -1,20 +0,0 @@ -#vertices of xfem subcells -0 0 -1 0 -1 1 -0 1 - -0.25 1 -1 0.25 - -0.625 0.625 -0.75 0.75 - -0.75 0.75 -0.625 1 - -0.75 0.75 -1 0.625 - -0 0 -0.625 0.625 diff --git a/examples/step-47/step-47.cc b/examples/step-47/step-47.cc deleted file mode 100644 index b5b98249fd..0000000000 --- a/examples/step-47/step-47.cc +++ /dev/null @@ -1,1128 +0,0 @@ -/* --------------------------------------------------------------------- - * - * Copyright (C) 2011 - 2018 by the deal.II authors - * - * This file is part of the deal.II library. - * - * The deal.II library is free software; you can use it, redistribute - * it, and/or modify it under the terms of the GNU Lesser General - * Public License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * The full text of the license can be found in the file LICENSE at - * the top level of the deal.II distribution. - * - * --------------------------------------------------------------------- - - * - * Author: Wolfgang Bangerth, University of Heidelberg, 2000 - */ - - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include - -#include -#include - - -#include - -#include - -#include - -namespace Step47 -{ - using namespace dealii; - - - - double sign (double d) - { - if (d > 0) - return 1; - else if (d < 0) - return -1; - else - return 0; - } - - - template - class LaplaceProblem - { - public: - LaplaceProblem (); - ~LaplaceProblem (); - - void run (); - - private: - 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); - - void setup_system (); - void assemble_system (); - void solve (); - void refine_grid (); - void output_results (const unsigned int cycle) const; - void compute_error () const; - - Triangulation triangulation; - - hp::DoFHandler dof_handler; - hp::FECollection fe_collection; - - ConstraintMatrix constraints; - - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; - - Vector solution; - Vector system_rhs; - }; - - - - - template - class Coefficient : public Function - { - public: - Coefficient () : Function() {} - - virtual double value (const Point &p, - const unsigned int component = 0) const override; - - virtual void value_list (const std::vector > &points, - std::vector &values, - const unsigned int component = 0) const override; - }; - - - - template - double Coefficient::value (const Point &p, - const unsigned int) const - { - if (p.square() < 0.5*0.5) - return 20; - else - return 1; - } - - - - template - void Coefficient::value_list (const std::vector > &points, - std::vector &values, - const unsigned int component) const - { - const unsigned int n_points = points.size(); - - Assert (values.size() == n_points, - ExcDimensionMismatch (values.size(), n_points)); - - (void) component; - Assert(component == 0, ExcIndexRange(component, 0, 1)); - - for (unsigned int i=0; i - double exact_solution (const Point &p) - { - const double r = p.norm(); - - return (r < 0.5 - ? - 1./20 * (-1./4*r*r + 61./16) - : - 1./4 * (1-r*r)); - } - - - template - LaplaceProblem::LaplaceProblem () - : - dof_handler (triangulation) - { - fe_collection.push_back (FESystem (FE_Q(1), 1, - FE_Nothing(), 1)); - fe_collection.push_back (FESystem (FE_Q(1), 1, - FE_Q(1), 1)); - } - - - - template - LaplaceProblem::~LaplaceProblem () - { - dof_handler.clear (); - } - - - - template - double - level_set (const Point &p) - { - return p.norm() - 0.5; - } - - - - template - Tensor<1,dim> - grad_level_set (const Point &p) - { - 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; - } - - - - template - void LaplaceProblem::setup_system () - { - for (typename hp::DoFHandler::cell_iterator cell - = dof_handler.begin_active(); - cell != dof_handler.end(); ++cell) - if (interface_intersects_cell(cell) == false) - cell->set_active_fe_index(0); - else - cell->set_active_fe_index(1); - - dof_handler.distribute_dofs (fe_collection); - - solution.reinit (dof_handler.n_dofs()); - system_rhs.reinit (dof_handler.n_dofs()); - - - constraints.clear (); -//TODO: fix this, it currently crashes - // DoFTools::make_hanging_node_constraints (dof_handler, constraints); - -//TODO: component 1 must satisfy zero boundary conditions - constraints.close(); - - - DynamicSparsityPattern dsp(dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern (dof_handler, dsp); - - constraints.condense (dsp); - - sparsity_pattern.copy_from(dsp); - - system_matrix.reinit (sparsity_pattern); - } - - - template - void LaplaceProblem::assemble_system () - { - const QGauss quadrature_formula(3); - - - FEValues plain_fe_values (fe_collection[0], quadrature_formula, - update_values | update_gradients | - update_quadrature_points | update_JxW_values); - - const unsigned int n_q_points = quadrature_formula.size(); - - FullMatrix cell_matrix; - Vector cell_rhs; - - std::vector local_dof_indices; - - const Coefficient coefficient; - std::vector coefficient_values (n_q_points); - - typename hp::DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - - 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); - - cell_matrix = 0; - cell_rhs = 0; - - if (cell->active_fe_index() == 0) - { - plain_fe_values.reinit (cell); - - 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()); - Assert (interface_intersects_cell(cell) == true, ExcInternalError()); - - 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)); - - FEValues this_fe_values (fe_collection[1], - compute_quadrature(quadrature_formula, cell, - level_set_values).second, - update_values | update_gradients | - update_quadrature_points | update_JxW_values ); - - this_fe_values.reinit (cell); - - 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))) - - - std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(j).second))))* - 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))) - - - std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(i).second))))* - 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))) - - - std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(i).second))))* - 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))) - - - std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(j).second))))* - 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))) - - - std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(i).second))))* - 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); - constraints.distribute_local_to_global (cell_matrix, cell_rhs, - local_dof_indices, - system_matrix, system_rhs); - } - - - std::map boundary_values; - VectorTools::interpolate_boundary_values (dof_handler, - 0, - Functions::ZeroFunction(2), - boundary_values); - MatrixTools::apply_boundary_values (boundary_values, - system_matrix, - solution, - 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 -// discontinuity 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 - 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; - - // 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 following 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 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); - - if (type == 1) - return std::pair >(1, plain_quadrature); - - 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); - - 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 - - // 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) = 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); - 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}} - }; - - 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); - } - } - - Assert (xfem_quadrature.size() == plain_quadrature.size() * 5, ExcInternalError()); - return std::pair >(2, xfem_quadrature); - } - - // 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(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); - 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++) - { - 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); - } - } - Assert (xfem_quadrature.size() == plain_quadrature.size() * 2, ExcInternalError()); - return std::pair >(3, xfem_quadrature); - } - - return std::pair >(0, plain_quadrature);; - - } - - 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. - - 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); - - const unsigned int n_q_points = plain_quadrature.size(); - - std::vector JxW(n_q_points); - - for ( unsigned int i = 0; i < n_q_points; i++) - { - 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; - } - - default: - Assert (false, ExcNotImplemented()); - } - - - Tensor<2,dim> jacobian; - - // 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); - } - - } - - - template - void LaplaceProblem::solve () - { - SolverControl solver_control (1000, 1e-12); - SolverCG<> solver (solver_control); - - PreconditionSSOR<> preconditioner; - preconditioner.initialize(system_matrix, 1.2); - - solver.solve (system_matrix, solution, system_rhs, - preconditioner); - - constraints.distribute (solution); - } - - - - template - void LaplaceProblem::refine_grid () - { - Vector estimated_error_per_cell (triangulation.n_active_cells()); - - KellyErrorEstimator::estimate (dof_handler, - QGauss(3), - typename FunctionMap::type(), - solution, - estimated_error_per_cell); - - GridRefinement::refine_and_coarsen_fixed_number (triangulation, - estimated_error_per_cell, - 0.3, 0.03); - - triangulation.execute_coarsening_and_refinement (); - } - - - - template - class Postprocessor : public DataPostprocessor - { - public: - virtual - void - evaluate_vector_field - (const dealii::DataPostprocessorInputs::Vector &inputs, - std::vector > &computed_quantities) const override; - - virtual std::vector get_names () const override; - - virtual - std::vector - get_data_component_interpretation () const override; - - virtual UpdateFlags get_needed_update_flags () const override; - }; - - - template - std::vector - Postprocessor::get_names() const - { - std::vector solution_names (1, "total_solution"); - solution_names.emplace_back("error"); - return solution_names; - } - - - template - std::vector - Postprocessor:: - get_data_component_interpretation () const - { - std::vector - interpretation (2, - DataComponentInterpretation::component_is_scalar); - return interpretation; - } - - - template - UpdateFlags - Postprocessor::get_needed_update_flags() const - { - return update_values | update_q_points; - } - - - template - void - Postprocessor:: - evaluate_vector_field - (const dealii::DataPostprocessorInputs::Vector &inputs, - std::vector > &computed_quantities) const - { - const unsigned int n_quadrature_points = inputs.solution_values.size(); - Assert (computed_quantities.size() == n_quadrature_points, - ExcInternalError()); - Assert (inputs.solution_values[0].size() == 2, - ExcInternalError()); - - for (unsigned int q=0; q - void LaplaceProblem::output_results (const unsigned int cycle) const - { - Assert (cycle < 10, ExcNotImplemented()); - - std::string filename = "solution-"; - filename += ('0' + cycle); - filename += ".vtk"; - - std::ofstream output (filename); - - Postprocessor postprocessor; - DataOut > data_out; - - data_out.attach_dof_handler (dof_handler); - data_out.add_data_vector (solution, "solution"); - data_out.add_data_vector (solution, postprocessor); - data_out.build_patches (5); - - data_out.write_vtk (output); - } - - - - template - void LaplaceProblem::compute_error () const - { - hp::QCollection q_collection; - q_collection.push_back (QGauss(2)); - q_collection.push_back (QIterated(QGauss<1>(2), 4)); - - hp::FEValues hp_fe_values (fe_collection, q_collection, - update_values | update_q_points | update_JxW_values); - - double l2_error_square = 0; - - std::vector > solution_values; - - typename hp::DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - - for (; cell!=endc; ++cell) - { - hp_fe_values.reinit (cell); - - const FEValues &fe_values = hp_fe_values.get_present_fe_values (); - - solution_values.resize (fe_values.n_quadrature_points, - Vector(2)); - fe_values.get_function_values (solution, - solution_values); - - for (unsigned int q=0; q laplace_problem_2d; - laplace_problem_2d.run (); - } - catch (std::exception &exc) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - - return 1; - } - catch (...) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } - - return 0; -}