Dear all,
I would be most grateful for your help on the following matter. I am
solving a combined stokes-laplace problem (have a vector valued solution).
I am trying to calculate the hdiv norm on the first two components of my
solution (i.e. the two components of the velocity), cellwise. The
implementation of my combined problem is from tutorial 46. The
integrate_difference is from tutorial 20 since I needed the mask for the
velocity components. However, I get an error saying 'no matching function
for the call to integrate difference'. Admittedly, the version of the
function i am using (i.e. with the component mask as the last included
argument) is not in vector_tools.h which is where all the other versions
of the integrate_differnce are I think. I include the file for your
consideration. I would be most grateful for your input. For succesful
implementations of this version of integrate_difference please refer to
tutorials 20, 55 and 56.
Kind regards,
Georgios
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
For more options, visit https://groups.google.com/d/optout.
/* ---------------------------------------------------------------------
*
* Copyright (C) 2011 - 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.
*
* ---------------------------------------------------------------------
*
* Author: Wolfgang Bangerth, Texas A&M University, 2011
*/
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/function.h>
#include <deal.II/base/utilities.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/constraint_matrix.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/hp/dof_handler.h>
#include <deal.II/hp/fe_collection.h>
#include <deal.II/hp/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/base/convergence_table.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include <deal.II/lac/block_vector.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/lac/block_sparse_matrix.h>
namespace Step46
{
using namespace dealii;
template <int dim>
class FluidStructureProblem
{
public:
FluidStructureProblem (const unsigned int stokes_degree,
const unsigned int laplace_degree);
void run ();
private:
enum
{
stokes_domain_id,
laplace_domain_id
};
static bool
cell_is_in_stokes_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell);
static bool
cell_is_in_laplace_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell);
void make_grid ();
void set_active_fe_indices ();
void setup_dofs ();
void assemble_system ();
void assemble_interface_term (const FEFaceValuesBase<dim> &laplace_fe_face_values,
const FEFaceValuesBase<dim> &stokes_fe_face_values,
std::vector<Tensor<1,dim> > &laplace_phi,
std::vector<Tensor<2,dim> > &laplace_grad_phi,
std::vector<SymmetricTensor<2,dim> > &stokes_symgrad_phi_u,
std::vector<double> &stokes_phi_p,
FullMatrix<double> &local_interface_matrix) const;
// void assemble_interface_term (const FEFaceValuesBase<dim> &elasticity_fe_face_values,
// const FEFaceValuesBase<dim> &stokes_fe_face_values,
// std::vector<Tensor<1,dim> > &elasticity_phi,
// std::vector<SymmetricTensor<2,dim> > &stokes_symgrad_phi_u,
// std::vector<double> &stokes_phi_p,
// FullMatrix<double> &local_interface_matrix) const;
void solve ();
void output_results (const unsigned int refinement_cycle) const;
void compute_errors (const unsigned int refinement_cycle);
void refine_mesh ();
const unsigned int stokes_degree;
const unsigned int laplace_degree;
Triangulation<dim> triangulation;
FESystem<dim> stokes_fe;
FESystem<dim> laplace_fe;
hp::FECollection<dim> fe_collection;
hp::DoFHandler<dim> dof_handler;
ConstraintMatrix constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
ConvergenceTable convergence_table;
Vector<double> solution;
Vector<double> system_rhs;
const double viscosity;
// const double lambda;
// const double mu;
};
template <int dim>
class StokesBoundaryValues : public Function<dim>
{
public:
StokesBoundaryValues () : Function<dim>(dim+1) {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
virtual void vector_value (const Point<dim> &p,
Vector<double> &value) const;
};
template <int dim>
double
StokesBoundaryValues<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
Assert (component < this->n_components,
ExcIndexRange (component, 0, this->n_components));
if (component == dim-1)
switch (dim)
{
case 2:
return 0; //std::sin(numbers::PI*p[0]);
case 3:
return 0; //std::sin(numbers::PI*p[0]) * std::sin(numbers::PI*p[1]);
default:
Assert (false, ExcNotImplemented());
}
return 0;
}
template <int dim>
void
StokesBoundaryValues<dim>::vector_value (const Point<dim> &p,
Vector<double> &values) const
{
for (unsigned int c=0; c<this->n_components; ++c)
values(c) = 0;// StokesBoundaryValues<dim>::value (p, c);
}
template <int dim>
class RightHandSide : public Function<dim>
{
public:
RightHandSide () : Function<dim>(dim+1) {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
virtual void vector_value (const Point<dim> &p,
Vector<double> &value) const;
};
template <int dim>
double
RightHandSide<dim>::value (const Point<dim> &/*p*/,
const unsigned int /*component*/) const
{
return 0;
}
template <int dim>
void
RightHandSide<dim>::vector_value (const Point<dim> &p,
Vector<double> &values) const
{
// values(0) = 1;
// values(1) = 1;
// values(2) = 1;
if (p[dim-1]>=0.5){
values(0) = - 400.0 * (pow(p[0],4) * (6*p[1] - 3) + pow(p[0],3) * (6 - 12 * p[1])
+ 3 * pow(p[0],2) * (4 * pow(p[1],3) - 6 * pow(p[1],2) + 4* p[1] -1)
- 6 * p[0] * p[1] * (2 * pow(p[1],2) - 3 * p[1] + 1)
+ p[1] * (2 * pow(p[1],2) - 3 * p[1] + 1))
+ 30 * pow(p[0] - 0.5,2)* pow(p[1]-0.5,2) - 3 * pow(1-p[0] , 2) * pow(p[1] - 0.5, 3);
values(1) = + 400.0 * (pow(p[1],4) * (6*p[0] - 3) + pow(p[1],3) * (6 - 12 * p[0])
+ 3 * pow(p[1],2) * (4 * pow(p[0],3) - 6 * pow(p[0],2) + 4* p[0] -1)
- 6 * p[1] * p[0] * (2 * pow(p[0],2) - 3 * p[0] + 1)
+ p[0] * (2 * pow(p[0],2) - 3 * p[0] + 1))
+ 20 * pow(p[0] - 0.5,3)* (p[1]-0.5) + 3 * pow(1-p[0] , 3) * pow(p[1] - 0.5, 2);
values(2) = 0.0;
}
else {
values(0) = - 400.0 * (pow(p[0],4) * (6*p[1] - 3) + pow(p[0],3) * (6 - 12 * p[1])
+ 3 * pow(p[0],2) * (4 * pow(p[1],3) - 6 * pow(p[1],2) + 4* p[1] -1)
- 6 * p[0] * p[1] * (2 * pow(p[1],2) - 3 * p[1] + 1)
+ p[1] * (2 * pow(p[1],2) - 3 * p[1] + 1));
//+ 30 * pow(p[0] - 0.5,2)* pow(p[1]-0.5,2) - 3 * pow(1-p[0] , 2) * pow(p[1] - 0.5, 3);
values(1) = + 400.0 * (pow(p[1],4) * (6*p[0] - 3) + pow(p[1],3) * (6 - 12 * p[0])
+ 3 * pow(p[1],2) * (4 * pow(p[0],3) - 6 * pow(p[0],2) + 4* p[0] -1)
- 6 * p[1] * p[0] * (2 * pow(p[0],2) - 3 * p[0] + 1)
+ p[0] * (2 * pow(p[0],2) - 3 * p[0] + 1));
//+ 20 * pow(p[0] - 0.5,3)* (p[1]-0.5) + 3 * pow(1-p[0] , 3) * pow(p[1] - 0.5, 2);
values(2) = 0.0;
}
}
template <int dim>
FluidStructureProblem<dim>::
FluidStructureProblem (const unsigned int stokes_degree,
const unsigned int laplace_degree)
:
stokes_degree (stokes_degree),
laplace_degree (laplace_degree),
triangulation (Triangulation<dim>::maximum_smoothing),
stokes_fe (FE_Q<dim>(stokes_degree+1), dim,
FE_Q<dim>(stokes_degree), 1),
laplace_fe (FE_Q<dim>(stokes_degree+1), dim,
FE_Nothing<dim>(), 1),
dof_handler (triangulation),
viscosity (1)//,
// lambda (1),
// mu (1)
{
fe_collection.push_back (stokes_fe);
fe_collection.push_back (laplace_fe);
}
template <int dim>
bool
FluidStructureProblem<dim>::
cell_is_in_stokes_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell)
{
return (cell->material_id() == stokes_domain_id);
}
template <int dim>
bool
FluidStructureProblem<dim>::
cell_is_in_laplace_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell)
{
return (cell->material_id() == laplace_domain_id);
}
template <int dim>
void
FluidStructureProblem<dim>::make_grid ()
{
GridGenerator::subdivided_hyper_cube (triangulation, 8, 0, 1);
for (typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active();
cell != triangulation.end(); ++cell)
for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
if (cell->face(f)->at_boundary()
&&
(cell->face(f)->center()[dim-1] >= 0.0))
cell->face(f)->set_all_boundary_ids(1);
for (typename Triangulation<dim>::active_cell_iterator
cell = dof_handler.begin_active();
cell != dof_handler.end(); ++cell)
if (cell->center()[dim-1] >= 0.5)
cell->set_material_id (stokes_domain_id);
else
cell->set_material_id (laplace_domain_id);
}
template <int dim>
void
FluidStructureProblem<dim>::set_active_fe_indices ()
{
for (typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active();
cell != dof_handler.end(); ++cell)
{
if (cell_is_in_stokes_domain(cell))
cell->set_active_fe_index (0);
else if (cell_is_in_laplace_domain(cell))
cell->set_active_fe_index (1);
else
Assert (false, ExcNotImplemented());
}
}
template <int dim>
void
FluidStructureProblem<dim>::setup_dofs ()
{
set_active_fe_indices ();
dof_handler.distribute_dofs (fe_collection);
{
constraints.clear ();
DoFTools::make_hanging_node_constraints (dof_handler,
constraints);
const FEValuesExtractors::Vector velocities(0);
VectorTools::interpolate_boundary_values (dof_handler,
1,
StokesBoundaryValues<dim>(),
constraints,
fe_collection.component_mask(velocities));
// const FEValuesExtractors::Vector displacements(dim+1);
// VectorTools::interpolate_boundary_values (dof_handler,
// 0,
// ZeroFunction<dim>(dim+1+),
// constraints,
// fe_collection.component_mask(displacements));
}
// {
// std::vector<types::global_dof_index> local_face_dof_indices (stokes_fe.dofs_per_face);
// for (typename hp::DoFHandler<dim>::active_cell_iterator
// cell = dof_handler.begin_active();
// cell != dof_handler.end(); ++cell)
// if (cell_is_in_fluid_domain (cell))
// for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
// if (!cell->at_boundary(f))
// {
// bool face_is_on_interface = false;
// if ((cell->neighbor(f)->has_children() == false)
// &&
// (cell_is_in_solid_domain (cell->neighbor(f))))
// face_is_on_interface = true;
// else if (cell->neighbor(f)->has_children() == true)
// {
// for (unsigned int sf=0; sf<cell->face(f)->n_children(); ++sf)
// if (cell_is_in_solid_domain (cell->neighbor_child_on_subface
// (f, sf)))
// {
// face_is_on_interface = true;
// break;
// }
// }
// if (face_is_on_interface)
// {
// cell->face(f)->get_dof_indices (local_face_dof_indices, 0);
// for (unsigned int i=0; i<local_face_dof_indices.size(); ++i)
// if (stokes_fe.face_system_to_component_index(i).first < dim)
// constraints.add_line (local_face_dof_indices[i]);
// }
// }
// }
constraints.close ();
std::cout << " Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl
<< " Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< std::endl;
{
DynamicSparsityPattern dsp (dof_handler.n_dofs(),
dof_handler.n_dofs());
Table<2,DoFTools::Coupling> cell_coupling (fe_collection.n_components(),
fe_collection.n_components());
Table<2,DoFTools::Coupling> face_coupling (fe_collection.n_components(),
fe_collection.n_components());
for (unsigned int c=0; c<fe_collection.n_components(); ++c)
for (unsigned int d=0; d<fe_collection.n_components(); ++d)
{
if (((c<dim+1) && (d<dim+1)
&& !((c==dim) && (d==dim))))
cell_coupling[c][d] = DoFTools::always;
if (((c<dim+1) && (d<dim+1)
&& !((c==dim) && (d==dim))))
face_coupling[c][d] = DoFTools::always;
}
// DoFTools::make_flux_sparsity_pattern (dof_handler, dsp);
DoFTools::make_flux_sparsity_pattern (dof_handler, dsp,
cell_coupling, face_coupling);
constraints.condense (dsp);
sparsity_pattern.copy_from (dsp);
}
system_matrix.reinit (sparsity_pattern);
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
}
template <int dim>
void FluidStructureProblem<dim>::assemble_system ()
{
system_matrix=0;
system_rhs=0;
const QGauss<dim> stokes_quadrature(stokes_degree+2);
const QGauss<dim> laplace_quadrature(laplace_degree+2);
hp::QCollection<dim> q_collection;
q_collection.push_back (stokes_quadrature);
q_collection.push_back (laplace_quadrature);
hp::FEValues<dim> hp_fe_values (fe_collection, q_collection,
update_values |
update_second_derivatives |
update_quadrature_points |
update_JxW_values |
update_gradients);
const QGauss<dim-1> common_face_quadrature(std::max (stokes_degree+2,
laplace_degree+2));
FEFaceValues<dim> stokes_fe_face_values (stokes_fe,
common_face_quadrature,
update_JxW_values |
update_values |
update_normal_vectors |
update_gradients);
FEFaceValues<dim> laplace_fe_face_values (laplace_fe,
common_face_quadrature,
update_values |
update_normal_vectors |
update_gradients);
FESubfaceValues<dim> stokes_fe_subface_values (stokes_fe,
common_face_quadrature,
update_JxW_values |
update_normal_vectors |
update_gradients);
FESubfaceValues<dim> laplace_fe_subface_values (laplace_fe,
common_face_quadrature,
update_values |
update_normal_vectors |
update_gradients);
const unsigned int stokes_dofs_per_cell = stokes_fe.dofs_per_cell;
const unsigned int laplace_dofs_per_cell = laplace_fe.dofs_per_cell;
const unsigned int stokes_n_q_points = stokes_quadrature.size();
const unsigned int laplace_n_q_points = laplace_quadrature.size();
FullMatrix<double> local_matrix;
FullMatrix<double> local_interface_matrix (laplace_dofs_per_cell,
stokes_dofs_per_cell);
Vector<double> local_rhs;
std::vector<types::global_dof_index> local_dof_indices;
std::vector<types::global_dof_index> neighbor_dof_indices (stokes_dofs_per_cell);
const RightHandSide<dim> right_hand_side;
std::vector<Vector<double>> rhs_values (stokes_n_q_points,
Vector<double>(dim+1)); // I put this here
std::vector<Vector<double>> Laplacian_values (stokes_n_q_points,
Vector<double>(dim+1)); // I put this here
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
// const FEValuesExtractors::Vector displacements (dim+1);
std::vector<SymmetricTensor<2,dim> > stokes_symgrad_phi_u (stokes_dofs_per_cell);
std::vector<Tensor<2,dim> > stokes_grad_phi_u (laplace_dofs_per_cell);
std::vector<double> stokes_div_phi_u (stokes_dofs_per_cell);
std::vector<double> stokes_phi_p (stokes_dofs_per_cell);
std::vector<Tensor<2,dim> > laplace_grad_phi (laplace_dofs_per_cell);
// std::vector<double> laplace_div_phi (laplace_dofs_per_cell);
std::vector<Tensor<1,dim> > laplace_phi (laplace_dofs_per_cell);
typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
hp_fe_values.reinit (cell);
const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values();
local_matrix.reinit (cell->get_fe().dofs_per_cell,
cell->get_fe().dofs_per_cell);
local_rhs.reinit (cell->get_fe().dofs_per_cell);
if (cell_is_in_stokes_domain (cell))
{
const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
Assert (dofs_per_cell == stokes_dofs_per_cell,
ExcInternalError());
right_hand_side.vector_value_list(fe_values.get_quadrature_points(), rhs_values);
for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
{
for (unsigned int k=0; k<dofs_per_cell; ++k)
{
stokes_symgrad_phi_u[k] = fe_values[velocities].symmetric_gradient (k, q);
stokes_div_phi_u[k] = fe_values[velocities].divergence (k, q);
stokes_phi_p[k] = fe_values[pressure].value (k, q);
// Laplacian_values[k] = fe_values[velocities].get
}
for (unsigned int i=0; i<dofs_per_cell; ++i){
for (unsigned int j=0; j<dofs_per_cell; ++j){
local_matrix(i,j) += (2 * viscosity * stokes_symgrad_phi_u[i] * stokes_symgrad_phi_u[j]
- stokes_div_phi_u[i] * stokes_phi_p[j]
- stokes_phi_p[i] * stokes_div_phi_u[j])
* fe_values.JxW(q);
}
const unsigned int component_i =
stokes_fe.system_to_component_index(i).first;//fe_values.shape_value(i,q)
if (component_i<dim){
double r = rhs_values[q](component_i);
local_rhs(i) += fe_values[velocities].value(i,q)[component_i] *
rhs_values[q](component_i) *
fe_values.JxW(q);
}
}
}
}
else
{
const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
Assert (dofs_per_cell == laplace_dofs_per_cell,
ExcInternalError());
right_hand_side.vector_value_list(fe_values.get_quadrature_points(), rhs_values);
for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
{
for (unsigned int k=0; k<dofs_per_cell; ++k)
{
laplace_grad_phi[k] = fe_values[velocities].gradient (k, q);
// elasticity_div_phi[k] = fe_values[displacements].divergence (k, q);
}
for (unsigned int i=0; i<dofs_per_cell; ++i){
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
local_matrix(i,j)
+= (
scalar_product(laplace_grad_phi[i], laplace_grad_phi[j])
)
*
fe_values.JxW(q);
// local_matrix(i,j)
// += (
// lambda *
// elasticity_div_phi[i] * elasticity_div_phi[j]
// +
// mu *
// scalar_product(elasticity_grad_phi[i], elasticity_grad_phi[j])
// +
// mu *
// scalar_product(elasticity_grad_phi[i], transpose(elasticity_grad_phi[j]))
// )
// *
// fe_values.JxW(q);
}
const unsigned int component_i =
laplace_fe.system_to_component_index(i).first;
local_rhs(i) += fe_values[velocities].value(i,q)[component_i] *
rhs_values[q](component_i) *
fe_values.JxW(q);
}
}
}
local_dof_indices.resize (cell->get_fe().dofs_per_cell);
cell->get_dof_indices (local_dof_indices);
constraints.distribute_local_to_global (local_matrix, local_rhs,
local_dof_indices,
system_matrix, system_rhs);
if (cell_is_in_laplace_domain (cell))
for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
if (cell->at_boundary(f) == false)
{
if ((cell->neighbor(f)->level() == cell->level())
&&
(cell->neighbor(f)->has_children() == false)
&&
cell_is_in_stokes_domain (cell->neighbor(f)))
{
laplace_fe_face_values.reinit (cell, f);
stokes_fe_face_values.reinit (cell->neighbor(f),
cell->neighbor_of_neighbor(f));
assemble_interface_term (laplace_fe_face_values, stokes_fe_face_values,
laplace_phi, laplace_grad_phi, stokes_symgrad_phi_u, stokes_phi_p,
local_interface_matrix);
cell->neighbor(f)->get_dof_indices (neighbor_dof_indices);
constraints.distribute_local_to_global(local_interface_matrix,
local_dof_indices,
neighbor_dof_indices,
system_matrix);
}
else if ((cell->neighbor(f)->level() == cell->level())
&&
(cell->neighbor(f)->has_children() == true))
{
for (unsigned int subface=0;
subface<cell->face(f)->n_children();
++subface)
if (cell_is_in_stokes_domain (cell->neighbor_child_on_subface
(f, subface)))
{
laplace_fe_subface_values.reinit (cell,
f,
subface);
stokes_fe_face_values.reinit (cell->neighbor_child_on_subface (f, subface),
cell->neighbor_of_neighbor(f));
assemble_interface_term (laplace_fe_subface_values,
stokes_fe_face_values,
laplace_phi, laplace_grad_phi,
stokes_symgrad_phi_u, stokes_phi_p,
local_interface_matrix);
cell->neighbor_child_on_subface (f, subface)
->get_dof_indices (neighbor_dof_indices);
constraints.distribute_local_to_global(local_interface_matrix,
local_dof_indices,
neighbor_dof_indices,
system_matrix);
}
}
else if (cell->neighbor_is_coarser(f)
&&
cell_is_in_stokes_domain(cell->neighbor(f)))
{
laplace_fe_face_values.reinit (cell, f);
stokes_fe_subface_values.reinit (cell->neighbor(f),
cell->neighbor_of_coarser_neighbor(f).first,
cell->neighbor_of_coarser_neighbor(f).second);
assemble_interface_term (laplace_fe_face_values,
stokes_fe_subface_values,
laplace_phi, laplace_grad_phi,
stokes_symgrad_phi_u, stokes_phi_p,
local_interface_matrix);
cell->neighbor(f)->get_dof_indices (neighbor_dof_indices);
constraints.distribute_local_to_global(local_interface_matrix,
local_dof_indices,
neighbor_dof_indices,
system_matrix);
}
}
}
}
template <int dim>
void
FluidStructureProblem<dim>::
assemble_interface_term (const FEFaceValuesBase<dim> &laplace_fe_face_values,
const FEFaceValuesBase<dim> &stokes_fe_face_values,
std::vector<Tensor<1,dim> > &laplace_phi,
std::vector<Tensor<2,dim> > &laplace_grad_phi,
std::vector<SymmetricTensor<2,dim> > &stokes_symgrad_phi_u,
std::vector<double> &stokes_phi_p,
FullMatrix<double> &local_interface_matrix) const
{
Assert (stokes_fe_face_values.n_quadrature_points ==
laplace_fe_face_values.n_quadrature_points,
ExcInternalError());
const unsigned int n_face_quadrature_points
= laplace_fe_face_values.n_quadrature_points;
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
//const FEValuesExtractors::Vector displacements (dim+1);
local_interface_matrix = 0;
for (unsigned int q=0; q<n_face_quadrature_points; ++q)
{
const Tensor<1,dim> normal_vector_stokes = stokes_fe_face_values.normal_vector(q);
const Tensor<1,dim> normal_vector_laplace = laplace_fe_face_values.normal_vector(q);
for (unsigned int k=0; k<stokes_fe_face_values.dofs_per_cell; ++k){
stokes_symgrad_phi_u[k] = stokes_fe_face_values[velocities].symmetric_gradient (k, q);
// stokes_grad_phi_u[k] = stokes_fe_face_values[velocities].gradient (k, q);
// stokes_phi_p[k] = stokes_fe_face_values[pressure].value(k,q);
}
for (unsigned int k=0; k<laplace_fe_face_values.dofs_per_cell; ++k){
laplace_phi[k] = laplace_fe_face_values[velocities].value (k,q);
// laplace_grad_phi[k] = laplace_fe_face_values[velocities].gradient(k,q);
}
for (unsigned int i=0; i<laplace_fe_face_values.dofs_per_cell; ++i)
for (unsigned int j=0; j<stokes_fe_face_values.dofs_per_cell; ++j)
local_interface_matrix(i,j) += 0;
// -((
// (stokes_symgrad_phi_u[j] *
// normal_vector_stokes)
// -
// stokes_phi_p[j] *
// normal_vector_stokes
// +
// laplace_grad_phi[j] *
// normal_vector_laplace) *
// laplace_phi[i] *
// stokes_fe_face_values.JxW(q));
}
}
template <int dim>
void
FluidStructureProblem<dim>::solve ()
{
SparseDirectUMFPACK direct_solver;
direct_solver.initialize (system_matrix);
direct_solver.vmult (solution, system_rhs);
constraints.distribute (solution);
}
template <int dim>
void
FluidStructureProblem<dim>::
output_results (const unsigned int refinement_cycle) const
{
std::vector<std::string> solution_names (dim, "velocity");
solution_names.push_back ("pressure");
// for (unsigned int d=0; d<dim; ++d)
// solution_names.push_back ("displacement");
std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_component_interpretation
(dim, DataComponentInterpretation::component_is_part_of_vector);
data_component_interpretation
.push_back (DataComponentInterpretation::component_is_scalar);
for (unsigned int d=0; d<dim; ++d)
data_component_interpretation
.push_back (DataComponentInterpretation::component_is_part_of_vector);
DataOut<dim,hp::DoFHandler<dim> > data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, solution_names,
DataOut<dim,hp::DoFHandler<dim> >::type_dof_data,
data_component_interpretation);
data_out.build_patches ();
std::ostringstream filename;
filename << "/Users/gs1511/dealii_out/solutiongs-"
<< Utilities::int_to_string (refinement_cycle, 2)
<< ".vtk";
std::ofstream output (filename.str().c_str());
data_out.write_vtk (output);
}
template <int dim>
void
FluidStructureProblem<dim>::refine_mesh ()
{
Vector<float>
stokes_estimated_error_per_cell (triangulation.n_active_cells());
Vector<float>
laplace_estimated_error_per_cell (triangulation.n_active_cells());
Vector<float>
model_error_per_cell (triangulation.n_active_cells());
Vector<float>
discretization_error_per_cell (triangulation.n_active_cells());
const QGauss<dim-1> stokes_face_quadrature(stokes_degree+2);
const QGauss<dim-1> laplace_face_quadrature(laplace_degree+2);
hp::QCollection<dim-1> face_q_collection;
face_q_collection.push_back (stokes_face_quadrature);
face_q_collection.push_back (laplace_face_quadrature);
const FEValuesExtractors::Vector velocities(0);
const FEValuesExtractors::Scalar pressure (dim);
KellyErrorEstimator<dim>::estimate (dof_handler,
face_q_collection,
typename FunctionMap<dim>::type(),
solution,
stokes_estimated_error_per_cell,
fe_collection.component_mask(velocities));
//const FEValuesExtractors::Vector velocities(0);
KellyErrorEstimator<dim>::estimate (dof_handler,
face_q_collection,
typename FunctionMap<dim>::type(),
solution,
laplace_estimated_error_per_cell,
fe_collection.component_mask(velocities));
stokes_estimated_error_per_cell
*= 4. / stokes_estimated_error_per_cell.l2_norm();
laplace_estimated_error_per_cell
*= 4. / laplace_estimated_error_per_cell.l2_norm();
Vector<float>
estimated_error_per_cell (triangulation.n_active_cells());
estimated_error_per_cell += stokes_estimated_error_per_cell;
estimated_error_per_cell += laplace_estimated_error_per_cell;
// // Make the terms for the laplace a-posteriori indicator (see pg 243 of Verfurth).
//
system_matrix=0;
system_rhs=0;
const QGauss<dim> stokes_quadrature(stokes_degree+2);
const QGauss<dim> laplace_quadrature(laplace_degree+2);
hp::QCollection<dim> q_collection;
q_collection.push_back (stokes_quadrature);
q_collection.push_back (laplace_quadrature);
hp::FEValues<dim> hp_fe_values (fe_collection, q_collection,
update_values |
update_quadrature_points |
update_JxW_values |
update_gradients);
const QGauss<dim-1> common_face_quadrature(std::max (stokes_degree+2,
laplace_degree+2));
// cell face calculations
FEFaceValues<dim> stokes_fe_face_values (stokes_fe,
common_face_quadrature,
update_JxW_values |
update_values |
update_normal_vectors |
update_gradients);
FEFaceValues<dim> laplace_fe_face_values (laplace_fe,
common_face_quadrature,
update_values |
update_normal_vectors |
update_gradients);
FESubfaceValues<dim> stokes_fe_subface_values (stokes_fe,
common_face_quadrature,
update_JxW_values |
update_normal_vectors |
update_gradients);
FESubfaceValues<dim> laplace_fe_subface_values (laplace_fe,
common_face_quadrature,
update_values |
update_normal_vectors |
update_gradients);
// neighbor of the cell face calculations
FEFaceValues<dim> stokes_fe_face_values_neighbor (stokes_fe,
common_face_quadrature,
update_JxW_values |
update_values |
update_normal_vectors |
update_gradients);
FEFaceValues<dim> laplace_fe_face_values_neighbor (laplace_fe,
common_face_quadrature,
update_values |
update_normal_vectors |
update_gradients);
FESubfaceValues<dim> stokes_fe_subface_values_neighbor (stokes_fe,
common_face_quadrature,
update_JxW_values |
update_normal_vectors |
update_gradients);
FESubfaceValues<dim> laplace_fe_subface_values_neighbor (laplace_fe,
common_face_quadrature,
update_values |
update_normal_vectors |
update_gradients);
const unsigned int stokes_dofs_per_cell = stokes_fe.dofs_per_cell;
const unsigned int laplace_dofs_per_cell = laplace_fe.dofs_per_cell;
const unsigned int stokes_n_q_points = stokes_quadrature.size();
const unsigned int laplace_n_q_points = laplace_quadrature.size();
const RightHandSide<dim> right_hand_side;
std::vector<Vector<double>> rhs_values (stokes_n_q_points,
Vector<double>(dim+1)); // I put this here
double cell_contribution;
double face_contribution;
std::vector<Tensor<2, dim> > cell_velocity_grads(stokes_n_q_points);
std::vector<Tensor<2, dim> > cell_neighbor_velocity_grads(stokes_n_q_points);
std::vector<Tensor<1,dim>> cell_laplacians(stokes_n_q_points);
std::vector<Tensor<1,dim>> cell_pressure_gradients(stokes_n_q_points);
std::vector<double> cell_values_pressure(stokes_n_q_points);
std::vector<double> cell_values_divergence(stokes_n_q_points);
std::vector<double> cell_neighbor_values_pressure(stokes_n_q_points);
// (1) : first term is ||f- Δu_T ||_{L^2}
for (typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active();
cell != dof_handler.end(); ++cell){
cell_contribution = 0;
hp_fe_values.reinit (cell);
const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values();
fe_values[pressure].get_function_gradients(solution, cell_pressure_gradients);
fe_values[velocities].get_function_laplacians(solution, cell_laplacians);
fe_values[velocities].get_function_divergences(solution, cell_values_divergence);
right_hand_side.vector_value_list(fe_values.get_quadrature_points(), rhs_values);
// for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q){
// for (unsigned int d = 0; d<dim; ++d){
// cell_contribution += pow(cell->diameter(),2) * pow(rhs_values[q][d] + cell_laplacians[q][d]-cell_pressure_gradients[q][d],2) * fe_values.JxW(q);
// }
// }
for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q){
cell_contribution += pow(cell_values_divergence[q] , 2) * fe_values.JxW(q);
}
// face stuff here
// for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no){
// face_contribution = 0;
// if (!( (cell_is_in_stokes_domain(cell) && cell_is_in_stokes_domain(cell->neighbor(face_no))) ||
// (cell_is_in_laplace_domain(cell) && cell_is_in_laplace_domain(cell->neighbor(face_no)))))
// {
// face_contribution = 0;
// }
//
//
// if (cell->face(face_no)->at_boundary())
// {
// face_integrals[cell->face(face_no)] = 0; // cell_contribution = 0;
// continue;
// }
// if ((cell->neighbor(face_no)->has_children() == false) &&
// (cell->neighbor(face_no)->level() == cell->level()) &&
// (cell->neighbor(face_no)->index() < cell->index()))
// continue;
// if (cell->at_boundary(face_no) == false)
// if (cell->neighbor(face_no)->level() < cell->level())
// continue;
//
//
// if (cell->face(face_no)->has_children() == false){
//// Integration over regular face (same level of refinement as the neighbour
// Assert (cell->neighbor(face_no).state() == IteratorState::valid, ExcInternalError());
//
// const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
// const active_cell_iterator neighbor = cell->neighbor(face_no);
//
// if (cell_is_in_stokes_domain(cell) && cell_is_in_stokes_domain(cell->neighbor(face_no))){
//
// stokes_fe_face_values.reinit(cell,face_no);
// stokes_fe_face_values[velocities].get_function_gradients(solution, cell_velocity_grads);
// stokes_fe_face_values[pressure].get_function_gradients(solution, cell_values_pressure);
//
// stokes_fe_face_values.reinit(cell->neighbor(face_no) , cell->neighbor_of_neighbor(face_no));
// stokes_fe_face_values[velocities].get_function_gradients(solution, cell_neighbor_velocity_grads);
// stokes_fe_face_values[pressure].get_function_gradients(solution, cell_neighbor_values_pressure);
//
// for (unsigned int q; q<stokes_fe_face_values.n_quadrature_points; ++q){
// const Tensor<1,dim> normal_vector_stokes = stokes_fe_face_values.normal_vector(q);
//
// face_contribution += pow(cell->diameter(),-0.5) * pow(((-cell_values_pressure[q] * normal_vector_stokes + cell_velocity_grads[q] * normal_vector_stokes)-(-cell_neighbor_values_pressure[q] * normal_vector_stokes + cell_neighbor_velocity_grads[q] * normal_vector_stokes)),2)* stokes_fe_face_values.JxW(q);
//
// }
// }
// else if (cell_is_in_laplace_domain(cell) && cell_is_in_laplace_domain(cell->neighbor(face_no))){
// laplace_fe_face_values.reinit(cell,face_no);
// laplace_fe_face_values[velocities].get_function_gradients(solution, cell_velocity_grads);
//
// laplace_fe_face_values.reinit(cell->neighbor(face_no) , cell->neighbor_of_neighbor(face_no));
// laplace_fe_face_values[velocities].get_function_gradients(solution, cell_neighbor_velocity_grads);
//
// for (unsigned int q; q<laplace_fe_face_values.n_quadrature_points; ++q){
// const Tensor<1,dim> normal_vector_laplace = laplace_fe_face_values.normal_vector(q);
//
// const Tensor<1,dim> normal_vector_laplace = laplace_fe_face_values.normal_vector(q);
// face_contribution += pow(cell->diameter(),-0.5) * pow(((cell_velocity_grads[q] * normal_vector_laplace)-(cell_neighbor_velocity_grads[q] * normal_vector_laplace)),2)* laplace_fe_face_values.JxW(q);
//
// }
// }
//
//
// }
//
// else
//// integrate_over_irregular_face (cell, face_no,
//// scratch_data.primal_solution,
//// scratch_data.dual_weights,
//// scratch_data.face_data,
//// face_integrals);
// }
// cell_contribution += face_contribution;
estimated_error_per_cell(cell->active_cell_index()) += (float) cell_contribution;
}
GridRefinement::refine_and_coarsen_fixed_number (triangulation,
estimated_error_per_cell,
0.3, 0.0);
triangulation.execute_coarsening_and_refinement ();
}
template <int dim>
void FluidStructureProblem<dim>::run ()
{
make_grid ();
for (unsigned int refinement_cycle = 0; refinement_cycle<8;
++refinement_cycle)
{
std::cout << "Refinement cycle " << refinement_cycle << std::endl;
if (refinement_cycle > 0)
refine_mesh ();
setup_dofs ();
std::cout << " Assembling..." << std::endl;
assemble_system ();
std::cout << " Solving..." << std::endl;
solve ();
std::cout << " Writing output..." << std::endl;
output_results (refinement_cycle);
compute_errors (refinement_cycle);
std::cout << std::endl;
}
}
template <int dim>
void FluidStructureProblem<dim>::compute_errors(const unsigned int refinement_cycle)
{
const ComponentSelectFunction<dim> pressure_mask (dim, dim+1);
const ComponentSelectFunction<dim> velocity_mask(std::make_pair(0, dim), dim+1);
Vector<double> cellwise_errors (triangulation.n_active_cells());
QTrapez<1> q_trapez;
QIterated<dim> quadrature (q_trapez, stokes_degree+2);
// VectorTools::integrate_difference (dof_handler, solution, exact_solution,
// cellwise_errors, quadrature,
// VectorTools::L2_norm,
// &pressure_mask);
// const double p_l2_error = VectorTools::compute_global_error(triangulation,
// cellwise_errors,
// VectorTools::L2_norm);
//
// VectorTools::integrate_difference (dof_handler, solution, exact_solution,
// cellwise_errors, quadrature,
// VectorTools::L2_norm,
// &velocity_mask);
// const double u_l2_error = VectorTools::compute_global_error(triangulation,
// cellwise_errors,
// VectorTools::L2_norm);
//
// VectorTools::integrate_difference (dof_handler, solution, exact_solution,
// cellwise_errors, quadrature,
// VectorTools::H1_norm,
// &velocity_mask);
// const double u_h1_error = VectorTools::compute_global_error(triangulation,
// cellwise_errors,
// VectorTools::H1_norm);
VectorTools::integrate_difference (dof_handler, solution, ZeroFunction<dim>(),
cellwise_errors, quadrature,
VectorTools::L2_norm,
&velocity_mask);
const double u_hdiv = VectorTools::compute_global_error(triangulation,
cellwise_errors,
VectorTools::Hdiv_seminorm);
const unsigned int n_active_cells=triangulation.n_active_cells();
const unsigned int n_dofs=dof_handler.n_dofs();
std::cout << "Cycle " << refinement_cycle << ':'
<< std::endl
<< " Number of active cells: "
<< n_active_cells
<< std::endl
<< " Number of degrees of freedom: "
<< n_dofs
<< std::endl;
std::cout << "Divergence of u_h: ||u_h||_Hdiv = " << u_hdiv
// << ", ||e_u||_L2 = " << u_l2_error
// << ", ||e_u||_H1 = " << u_h1_error
// << ", ||div(u_h)||_Hdiv = " << u_hdiv
<< std::endl;
convergence_table.add_value("refinement cycle", refinement_cycle);
convergence_table.add_value("cells", n_active_cells);
convergence_table.add_value("dofs", n_dofs);
// convergence_table.add_value("$||p-p_h||_{L^2}$", p_l2_error);
// convergence_table.add_value("$||u-u_h||_{L^2}$", u_l2_error);
// convergence_table.add_value("$||u-u_h||_{H^1}$", u_h1_error);
convergence_table.add_value("$||u_h||_{Hdiv}$", u_hdiv);
// convergence_table.set_scientific("$||p-p_h||_{L^2}$", true);
// convergence_table.set_scientific("$||u-u_h||_{L^2}$", true);
// convergence_table.set_scientific("$||u-u_h||_{H^1}$", true);
convergence_table.set_scientific("$||u_h||_{Hdiv}$", true);
std::cout << std::scientific<< std::endl;
convergence_table.write_text(std::cout<< std::scientific);
std::string error_filename = "/Users/gs1511/dealii_out/combined_gs_error-adaptive.tex";
std::ofstream error_table_file(error_filename.c_str());
convergence_table.write_tex(error_table_file);
std::string error_filename2 = "/Users/gs1511/dealii_out/combined_gs_error-adaptive.text";
std::ofstream error_table_file2(error_filename2.c_str());
convergence_table.write_text(error_table_file2);
}
}
int main ()
{
try
{
using namespace dealii;
using namespace Step46;
FluidStructureProblem<2> flow_problem(1, 1);
flow_problem.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;
}
