const FEValuesExtractors::Vector vec (first_vector_component);
- // coordinate directions of the face.
- const unsigned int
- face_coordinate_direction[GeometryInfo<dim>::faces_per_cell]
- = { 1, 1, 0, 0 };
-
- const double tol = 0.5 * cell->face (face)->diameter () / cell->get_fe ().degree;
- const std::vector<Point<dim> > &
- reference_quadrature_points = fe_face_values.get_quadrature_points ();
-
// Project the boundary function onto the shape functions for this edge
// and set up a linear system of equations to get the values for the DoFs
// associated with this edge.
for (unsigned int q_point = 0;
q_point < fe_face_values.n_quadrature_points; ++q_point)
{
- // Compute the tangent of the face
- // at the quadrature point.
- Point<dim> shifted_reference_point_1
- = reference_quadrature_points[q_point];
- Point<dim> shifted_reference_point_2
- = reference_quadrature_points[q_point];
-
- shifted_reference_point_1 (face_coordinate_direction[face])
- += tol;
- shifted_reference_point_2 (face_coordinate_direction[face])
- -= tol;
- Tensor<1,dim> tangential
- = (fe_face_values.get_mapping ()
- .transform_unit_to_real_cell (cell,
- shifted_reference_point_1)
- -
- fe_face_values.get_mapping ()
- .transform_unit_to_real_cell (cell,
- shifted_reference_point_2))
- / tol;
- tangential
- /= tangential.norm ();
-
// Compute the entires of the linear system
// Note the system is symmetric so we could only compute the lower/upper triangle.
//
//
// The RHS entries are:
// \int_{edge} (tangential* boundary_value) * (tangential * edge_shape_function_i) dS.
+ //
+ // In 2D, tangential*vector is equivalent to cross_product(normal, vector), so we use this instead.
+ // This avoids possible issues with the computation of the tangent.
+
+ // Store the normal at this quad point:
+ Point<dim> normal_at_q_point = fe_face_values.normal_vector(q_point);
for (unsigned int j = 0; j < associated_edge_dofs; ++j)
{
const unsigned int j_face_idx = associated_edge_dof_to_face_dof[j];
+ const unsigned int j_cell_idx = fe.face_to_cell_index (j_face_idx, face);
+
+ Tensor<1,dim> phi_j = fe_face_values[vec].value (j_cell_idx, q_point);
for (unsigned int i = 0; i < associated_edge_dofs; ++i)
{
const unsigned int i_face_idx = associated_edge_dof_to_face_dof[i];
+ const unsigned int i_cell_idx = fe.face_to_cell_index (i_face_idx, face);
+
+ Tensor<1,dim> phi_i = fe_face_values[vec].value (i_cell_idx, q_point);
+
+ // Using n cross phi
edge_matrix(i,j)
+= fe_face_values.JxW (q_point)
- * (fe_face_values[vec].value (fe.face_to_cell_index (i_face_idx, face), q_point) * tangential)
- * (fe_face_values[vec].value (fe.face_to_cell_index (j_face_idx, face), q_point) * tangential);
+ * ((phi_i[1]*normal_at_q_point(0) - phi_i[0]*normal_at_q_point(1))
+ * (phi_j[1]*normal_at_q_point(0) - phi_j[0]*normal_at_q_point(1)));
}
+ // Using n cross phi
edge_rhs(j)
+= fe_face_values.JxW (q_point)
- * (values[q_point] (first_vector_component) * tangential [0]
- + values[q_point] (first_vector_component + 1) * tangential [1])
- * (fe_face_values[vec].value (fe.face_to_cell_index (j_face_idx, face), q_point) * tangential);
+ * ((values[q_point] (first_vector_component+1) * normal_at_q_point (0)
+ - values[q_point] (first_vector_component) * normal_at_q_point (1))
+ * (phi_j[1]*normal_at_q_point(0) - phi_j[0]*normal_at_q_point(1)));
}
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2015 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.
+//
+// ---------------------------------------------------------------------
+//
+// By Ross Kynch
+//
+// Test to confirm that FE_Nedelec works on deformed elements in 2D when using
+// Dirichlet boundary conditions (n x E = f, type). This is handled by the
+// function project_boundary_values_curl_conforming_l2().
+//
+// This test solves the real valued curl-curl equation in 2D:
+//
+// curl(curl(E)) + E = Js
+//
+// where the solution is:
+// E = (0, x^2)
+//
+// so, Js = (0,-2) + E.
+//
+// i.e the solution should be "exact" at p=2.
+//
+// The domain is a distorted quad after 1 refinement. The distortion is
+// performed by GridTools::distort_random.
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/constraint_matrix.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_renumbering.h>
+
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+
+#include <fstream>
+#include <iostream>
+#include <sstream>
+
+
+using namespace dealii;
+
+namespace polytest
+{
+ template<int dim>
+ class SimplePolynomial : public Function<dim>
+ {
+ public:
+ SimplePolynomial() : Function<dim>() {}
+ void vector_value_list (const std::vector<Point<dim> > &points,
+ std::vector<Vector<double> > &values) const;
+
+ void rhs_value_list (const std::vector<Point<dim> > &points,
+ std::vector<Vector<double> > &values) const;
+
+ };
+ template<int dim>
+ void SimplePolynomial<dim>::vector_value_list (const std::vector<Point<dim> > &points,
+ std::vector<Vector<double> > &values) const
+ {
+ Assert (dim == 2, ExcNotImplemented());
+ Assert (values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
+ for (unsigned int i=0; i<points.size(); ++i)
+ {
+ const Point<dim> &p = points[i];
+ // non-zero curl-curl:
+ values[i][0] = 0.0;
+ values[i][1] = p[0]*p[0];
+
+ }
+
+ }
+
+ template<int dim>
+ void SimplePolynomial<dim>::rhs_value_list (const std::vector<Point<dim> > &points,
+ std::vector<Vector<double> > &values) const
+ {
+ Assert (dim == 2, ExcNotImplemented());
+ Assert (values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
+
+ for (unsigned int i=0; i<points.size(); ++i)
+ {
+ const Point<dim> &p = points[i];
+ // non-zero curl-curl:
+ values[i][0] = 0.0;
+ values[i][1] = -2.0 + p[0]*p[0];
+ }
+
+ }
+
+ template<int dim>
+ class polytest
+ {
+ public:
+ polytest(unsigned int degree);
+ ~polytest ();
+
+ void run();
+
+ private:
+ void setup_system();
+ void assemble_system();
+ void solve();
+
+ void output_error();
+
+ unsigned int p_order;
+ unsigned int quad_order;
+
+ Triangulation<dim> tria;
+ DoFHandler<dim> dof_handler;
+ FE_Nedelec<dim> fe;
+ ConstraintMatrix constraints;
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+ Vector<double> solution;
+ Vector<double> system_rhs;
+
+ };
+
+ template<int dim>
+ polytest<dim>::polytest (unsigned int degree)
+ :
+ p_order(degree),
+ quad_order(2*degree + 3),
+ dof_handler(tria),
+ fe (degree)
+ {
+ }
+ template <int dim>
+ polytest<dim>::~polytest ()
+ {
+ dof_handler.clear ();
+ }
+
+ template<int dim>
+ void polytest<dim>::setup_system()
+ {
+ dof_handler.distribute_dofs (fe);
+ solution.reinit (dof_handler.n_dofs());
+ system_rhs.reinit (dof_handler.n_dofs());
+
+
+ constraints.clear ();
+ SimplePolynomial<dim> boundary_function;
+
+ VectorTools::project_boundary_values_curl_conforming_l2 (dof_handler,
+ 0,
+ boundary_function,
+ 0,
+ constraints);
+ constraints.close ();
+ DynamicSparsityPattern c_sparsity(dof_handler.n_dofs());
+ DoFTools::make_sparsity_pattern(dof_handler,
+ c_sparsity,
+ constraints,false);
+
+ sparsity_pattern.copy_from(c_sparsity);
+ system_matrix.reinit (sparsity_pattern);
+ }
+
+ template<int dim>
+ void polytest<dim>::assemble_system()
+ {
+ const QGauss<dim> test_quad(quad_order);
+ FEValues<dim> fe_values_test (fe, test_quad,
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
+
+ const QGauss<dim> quadrature_formula(quad_order);
+ const unsigned int n_q_points = quadrature_formula.size();
+
+ const QGauss<dim-1> face_quadrature_formula(quad_order);
+ const unsigned int n_face_q_points = face_quadrature_formula.size();
+
+ FEValues<dim> fe_values (fe, quadrature_formula,
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
+
+ FEFaceValues<dim> fe_face_values(fe, face_quadrature_formula,
+ update_values | update_quadrature_points |
+ update_normal_vectors | update_JxW_values);
+
+ const FEValuesExtractors::Vector vec(0);
+
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+
+
+ FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs (dofs_per_cell);
+
+ std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+ SimplePolynomial<dim> right_hand_side;
+ std::vector<Vector<double> > rhs_value_list(n_q_points,
+ Vector<double>(fe.n_components()));
+
+ typename DoFHandler<dim>::active_cell_iterator cell, endc;
+ endc = dof_handler.end();
+ cell = dof_handler.begin_active();
+ for (; cell!=endc; ++cell)
+ {
+ fe_values_test.reinit(cell);
+ fe_values.reinit (cell);
+ cell_matrix=0;
+ cell_rhs=0;
+
+ right_hand_side.rhs_value_list(fe_values.get_quadrature_points(), rhs_value_list);
+ for (unsigned int q=0; q<n_q_points; ++q)
+ {
+ Tensor<1,dim> rhs_value;
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ rhs_value[d] = rhs_value_list[q](d);
+ }
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ cell_matrix(i,j)+= (
+ fe_values[vec].curl(i,q)*fe_values[vec].curl(j,q)
+ + fe_values[vec].value(i,q)*fe_values[vec].value(j,q)
+ )*fe_values.JxW(q);
+ }
+ cell_rhs(j) += rhs_value*fe_values[vec].value(j,q)*fe_values.JxW(q);
+ }
+ }
+ cell->get_dof_indices (local_dof_indices);
+ constraints.distribute_local_to_global(cell_matrix,
+ cell_rhs,
+ local_dof_indices,
+ system_matrix, system_rhs);
+ }
+ }
+
+ template<int dim>
+ void polytest<dim>::solve()
+ {
+ SparseDirectUMFPACK direct;
+ direct.initialize(system_matrix);
+
+ direct.vmult (solution, system_rhs);
+ constraints.distribute (solution);
+ }
+
+ template<int dim>
+ void polytest<dim>::output_error()
+ {
+ SimplePolynomial<dim> exact_solution;
+ Vector<double> diff_per_cell(tria.n_active_cells());
+ VectorTools::integrate_difference(dof_handler, solution, exact_solution,
+ diff_per_cell, QGauss<dim>(quad_order), VectorTools::L2_norm);
+ const double L2_error = diff_per_cell.l2_norm();
+
+ deallog << "p=" << p_order << " L2_error: " << L2_error << std::endl;
+ }
+
+ template<int dim>
+ void polytest<dim>::run()
+ {
+ GridGenerator::hyper_cube (tria, -1.2, 1.3);
+ // REFINE
+ tria.refine_global(1);
+
+ // DISTORT ALL:
+ GridTools::distort_random (0.2, tria, false);
+
+ setup_system();
+ assemble_system();
+ solve();
+ output_error();
+ }
+}
+int main ()
+{
+ const unsigned int dim(2);
+
+ std::ofstream logfile("output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ for (unsigned int p=0; p<3; ++p)
+ {
+ polytest::polytest<dim> poly(p);
+ poly.run();
+ }
+}
projects / dealii.git / commitdiff