+double sign (double d)
+{
+ if (d > 0)
+ return 1;
+ else if (d < 0)
+ return -1;
+ else
+ return 0;
+}
+
+
template <int dim>
class LaplaceProblem
{
void run ();
private:
+ double level_set (const Point<dim> &p) const;
+ Tensor<1,dim> grad_level_set (const Point<dim> &p) const;
+
bool interface_intersects_cell (const typename Triangulation<dim>::cell_iterator &cell) const;
- std::pair<unsigned int, Quadrature<dim> > compute_quadrature(const Quadrature<dim> &plain_quadrature, const typename hp::DoFHandler<dim>::active_cell_iterator &cell, const std::vector<double> &level_set_values);
- void append_quadrature(const Quadrature<dim> &plain_quadrature,
- const std::vector<Point<dim> > &v ,
- std::vector<Point<dim> > &xfem_points,
- std::vector<double> &xfem_weights);
+ std::pair<unsigned int, Quadrature<dim> > compute_quadrature(const Quadrature<dim> &plain_quadrature, const typename hp::DoFHandler<dim>::active_cell_iterator &cell, const std::vector<double> &level_set_values);
+ void append_quadrature(const Quadrature<dim> &plain_quadrature,
+ const std::vector<Point<dim> > &v ,
+ std::vector<Point<dim> > &xfem_points,
+ std::vector<double> &xfem_weights);
void setup_system ();
void assemble_system ();
+template <int dim>
+double
+LaplaceProblem<dim>::
+level_set (const Point<dim> &p) const
+{
+ return p.norm() - 0.5;
+}
+
+
+
+template <int dim>
+Tensor<1,dim>
+LaplaceProblem<dim>::
+grad_level_set (const Point<dim> &p) const
+{
+ return p / p.norm();
+}
+
+
+
template <int dim>
bool
LaplaceProblem<dim>::
interface_intersects_cell (const typename Triangulation<dim>::cell_iterator &cell) const
{
+ for (unsigned int v=0; v<GeometryInfo<dim>::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;
}
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
// constrain these DoFs
for (unsigned int v=0; v<GeometryInfo<dim>::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();
FEValues<dim> plain_fe_values (fe_collection[0], quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
+ FEValues<dim> 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();
cell = dof_handler.begin_active(),
endc = dof_handler.end();
- std::vector<double> 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<double> 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<unsigned int, Quadrature<dim> > type_and_quadrature = compute_quadrature(quadrature_formula, cell, level_set_values);
-
- std::cout << "type : " << type_and_quadrature.first << std::endl;
- std::vector<Point<dim> > points = type_and_quadrature.second.get_points();
- std::vector<double> 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; i<points.size(); i++)
- output << points[i] << std::endl;
- std::string filename2 = "weights.dat";
- std::ofstream output2 (filename2.c_str());
- output2 << "#xfem Weights" << std::endl;
- for (unsigned int i=0; i<weights.size(); i++)
- output2 << weights[i] << std::endl;
- assert(0);
-
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;
- 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_point<n_q_points; ++q_point)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- cell_matrix(i,j) += (coefficient_values[q_point] *
- plain_fe_values.shape_grad(i,q_point) *
- plain_fe_values.shape_grad(j,q_point) *
- plain_fe_values.JxW(q_point));
-
- cell_rhs(i) += (plain_fe_values.shape_value(i,q_point) *
- 1.0 *
- plain_fe_values.JxW(q_point));
- }
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ cell_matrix(i,j) += (coefficient_values[q_point] *
+ plain_fe_values.shape_grad(i,q_point) *
+ plain_fe_values.shape_grad(j,q_point) *
+ plain_fe_values.JxW(q_point));
+
+ cell_rhs(i) += (plain_fe_values.shape_value(i,q_point) *
+ 1.0 *
+ plain_fe_values.JxW(q_point));
+ }
+ }
+ else
+ {
+ Assert (cell->active_fe_index() == 1, ExcInternalError());
+
+ std::auto_ptr<FEValues<dim> > custom_fe_values;
+
+ if (interface_intersects_cell(cell) == false)
+ enriched_fe_values.reinit (cell);
+ else
+ {
+ std::vector<double> level_set_values (GeometryInfo<dim>::vertices_per_cell);
+ for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+ level_set_values[v] = level_set (cell->vertex(v));
+
+ custom_fe_values
+ .reset(new FEValues<dim> (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<dim> &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_point<this_fe_values.n_quadrature_points; ++q_point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ if (cell->get_fe().system_to_component_index(i).first == 0)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ if (cell->get_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; j<dofs_per_cell; ++j)
+ if (cell->get_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);
// are considered.
// Type 1: there is not cut. Type 2: a corner of the element is cut. Type 3: two corners are cut.
- template <int dim>
- std::pair<unsigned int, Quadrature<dim> >
- LaplaceProblem<dim>::compute_quadrature (const Quadrature<dim> &plain_quadrature,
- const typename hp::DoFHandler<dim>::active_cell_iterator &cell,
- const std::vector<double> &level_set_values )
+template <int dim>
+std::pair<unsigned int, Quadrature<dim> >
+LaplaceProblem<dim>::compute_quadrature (const Quadrature<dim> &plain_quadrature,
+ const typename hp::DoFHandler<dim>::active_cell_iterator &cell,
+ const std::vector<double> &level_set_values )
{
- unsigned int type = 0;
+ unsigned int type = 0;
- // find the type of cut
- int sign_ls[GeometryInfo<dim>::vertices_per_cell];
- for (unsigned int v=0; v<GeometryInfo<dim>::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<dim>::vertices_per_cell];
+ for (unsigned int v=0; v<GeometryInfo<dim>::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<dim> v0(0,0);
- Point<dim> v1(1,0);
- Point<dim> v2(0,1);
- Point<dim> v3(1,1);
+ Point<dim> v0(0,0);
+ Point<dim> v1(1,0);
+ Point<dim> v2(0,1);
+ Point<dim> v3(1,1);
- Point<dim> A(0,0);
- Point<dim> B(0,0);
- Point<dim> C(0,0);
- Point<dim> D(0,0);
- Point<dim> E(0,0);
- Point<dim> F(0,0);
+ Point<dim> A(0,0);
+ Point<dim> B(0,0);
+ Point<dim> C(0,0);
+ Point<dim> D(0,0);
+ Point<dim> E(0,0);
+ Point<dim> F(0,0);
- if (type == 1)
- return std::pair<unsigned int, Quadrature<dim> >(1, plain_quadrature);
+ if (type == 1)
+ return std::pair<unsigned int, Quadrature<dim> >(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<dim> xfem_quadrature(5*n_q_points);
+ Quadrature<dim> xfem_quadrature(5*n_q_points);
- std::vector<Point<dim> > v(GeometryInfo<dim>::vertices_per_cell);
+ std::vector<Point<dim> > v(GeometryInfo<dim>::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<dim> 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<Point<dim> > xfem_points;
- std::vector<double> 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<dim> 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<Point<dim> > xfem_points;
+ std::vector<double> 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<Point<dim> > 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<Point<dim> > 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<unsigned int, Quadrature<dim> >(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<unsigned int, Quadrature<dim> >(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<dim> xfem_quadrature(5*n_q_points);
+ // loop over all subelements for integration
+ // in type 2 there are 5 subelements
- std::vector<Point<dim> > v(GeometryInfo<dim>::vertices_per_cell);
+ Quadrature<dim> 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<Point<dim> > v(GeometryInfo<dim>::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<dim> 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<Point<dim> > xfem_points;
- std::vector<double> 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<dim> 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<Point<dim> > xfem_points;
+ std::vector<double> 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<Point<dim> > 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<Point<dim> > 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<unsigned int, Quadrature<dim> >(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<unsigned int, Quadrature<dim> >(3, xfem_quadrature);
+ }
- return std::pair<unsigned int, Quadrature<dim> >(0, plain_quadrature);;
+ return std::pair<unsigned int, Quadrature<dim> >(0, plain_quadrature);;
}
- template <int dim>
+template <int dim>
void LaplaceProblem<dim>::append_quadrature ( const Quadrature<dim> &plain_quadrature,
const std::vector<Point<dim> > &v,
std::vector<Point<dim> > &xfem_points,
std::vector<double> &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<dim>::vertices_per_cell;
+ unsigned int n_v = GeometryInfo<dim>::vertices_per_cell;
- std::vector<Point<dim> > q_points = plain_quadrature.get_points();
- std::vector<Point<dim> > q_transf(q_points.size());
- std::vector<double> W = plain_quadrature.get_weights();
- std::vector<double> phi(n_v);
- std::vector<Tensor<1,dim> > grad_phi(n_v);
+ std::vector<Point<dim> > q_points = plain_quadrature.get_points();
+ std::vector<Point<dim> > q_transf(q_points.size());
+ std::vector<double> W = plain_quadrature.get_weights();
+ std::vector<double> phi(n_v);
+ std::vector<Tensor<1,dim> > grad_phi(n_v);
- const unsigned int n_q_points = plain_quadrature.size();
+ const unsigned int n_q_points = plain_quadrature.size();
- std::vector<double> JxW(n_q_points);
+ std::vector<double> 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<dim; ++d)
- for (unsigned int e=0; e<dim; ++e)
- {
- for (unsigned int j = 0; j<GeometryInfo<dim>::vertices_per_cell; j++)
- {
+ // Calculate Jacobian of transformation
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int e=0; e<dim; ++e)
+ {
+ for (unsigned int j = 0; j<GeometryInfo<dim>::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<dim> q_prime;
- for (unsigned int d=0; d<dim; ++d)
- for (unsigned int j = 0; j<GeometryInfo<dim>::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<dim> q_prime;
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int j = 0; j<GeometryInfo<dim>::vertices_per_cell; j++)
+ q_prime[d] += v[j](d) * phi[j];
+ xfem_points.push_back(q_prime);
+ }
}
static const HyperBallBoundary<dim> boundary;
triangulation.set_boundary (0, boundary);
- triangulation.refine_global (1);
+ triangulation.refine_global (3);
}
else
refine_grid ();