Dear all,
I am working on a combined stokes laplace problem (I also attach the code).
Essentially I am solving stokes in one part of the domain and laplace in
another. To this extend, I want my velocity to be both continuous and and
have continuous first derivatives everywhere in my domain (including on
between interfaces). My implementation is similar to step 46. I am
choosing my spaces like I show below (in the stokes part and in the laplace
part of my domain respectively):
stokes_fe (FE_Q<dim>(stokes_degree+2), dim,
FE_Q<dim>(stokes_degree+1), 1),
laplace_fe (FE_Q<dim>(stokes_degree+2), dim,
FE_Nothing<dim>(), 1),
spaces with the intention of getting zero jump in the velocity gradient
between the interfaces, however, this does not happen. I have taken care
to choose my pressure in L^2_0 on both parts of the domain. Any suggestion
on how I can enforce continuity of velocity gradients across the interface?
I would be most grateful for your help.
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.
//
// stokes_laplace_schur.cpp
// stokes_cg
//
// Created by Georgios Sialounas on 19/08/2018.
// Copyright © 2018 Georgios Sialounas. All rights reserved.
//
/* ---------------------------------------------------------------------
*
* 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/block_vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.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/tria_boundary_lib.h>
//#include <deal.II/grid/grid_tools.h>
//#include <deal.II/grid/grid_refinement.h>
//#include <deal.II/grid/grid_out.h>
//
//#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
//#include <deal.II/dofs/dof_accessor.h>
//#include <deal.II/dofs/dof_tools.h>
//
//#include <deal.II/fe/fe_q.h>
//#include <deal.II/fe/fe_system.h>
//#include <deal.II/fe/fe_values.h>
//
//#include <deal.II/numerics/vector_tools.h>
//#include <deal.II/numerics/matrix_tools.h>
//#include <deal.II/numerics/data_out.h>
//#include <deal.II/numerics/error_estimator.h>
//#include <deal.II/base/convergence_table.h>
//// Then we need to include the header file for the sparse direct solver
//// UMFPACK:
//#include <deal.II/lac/sparse_direct.h>
//
//// This includes the library for the incomplete LU factorization that will be
//// used as a preconditioner in 3D:
#include <deal.II/lac/sparse_ilu.h>
///////////
//////////
//
//
//
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/block_sparse_matrix.h>
#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>
namespace Step46
{
using namespace dealii;
template <int dim>
struct InnerPreconditioner;
template <>
struct InnerPreconditioner<2>
{
typedef SparseDirectUMFPACK type;
};
template <>
struct InnerPreconditioner<3>
{
typedef SparseILU<double> type;
};
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;
//
// Vector<double> solution;
// Vector<double> system_rhs;
BlockSparsityPattern sparsity_pattern;
BlockSparseMatrix<double> system_matrix;
BlockVector<double> solution;
BlockVector<double> system_rhs;
std_cxx11::shared_ptr<typename InnerPreconditioner<dim>::type> A_preconditioner;
ConvergenceTable convergence_table;
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.0 ){
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.0)
- 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));
//+ (3* pow(p[0],2) - 3*p[0] + 0.5 ) * (pow(p[1], 3) - 1.5*pow(p[1],2) + 0.625*p[1]- 0.0625); // cubic pressure L20
//+ (3* pow(p[0],2) - 3*p[0] + 0.5 ) * (pow(p[1], 2) - p[1] + (double)(1.0/((double)6.0)));
//+ 30 * pow(p[0] - 0.5,2)* pow(p[1],2) - 3 * pow(1-p[0] , 2) * pow(p[1] - 0.5, 3);
//+ (3* pow(p[0],2) - 3*p[0] + 0.5 ) * (pow(p[1], 2) - p[1] + (double)(1.0/((double)6.0)));
//+ (3* pow(p[0],2) - 3*p[0] + 0.5 ) * (p[1] - 0.5)*(p[1] - 0.75) * (p[1] - 1.0);
//+ (3* pow(p[0],2) - 3*p[0] + 0.5 ) * (pow(p[1], 2) - p[1] + (double)(1.0/((double)6.0)));
//(pow(p[1],3) - 2.25*pow(p[1],2) +1.625*p[1] - 0.375);
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));
// + p[0] * (p[0] - 0.5) * (p[0] - 1) * (3*pow(p[1],2) -3*p[1]+0.625); // Cubic Pressure L20
// + p[0] * (p[0] - 0.5) * (p[0] - 1) * (2* p[1] - 1.0); //Quadratic pressure
// + 20 * pow(p[0] - 0.5,3)* p[1] + 3 * pow(1-p[0] , 3) * pow(p[1] - 0.5, 2);
//+ p[0] * (p[0] - 0.5) * (p[0] - 1) * (2* p[1] - 1.0);
// + p[0] * (p[0] - 0.5) * (p[0] - 1) * (3* pow(p[1], 2) - 4.5 * p[1] + 1.625);
// + p[0] * (p[0] - 0.5) * (p[0] - 1) * (2* p[1] - 1.0);
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],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>
class ExactSolution : public Function<dim>
{
public:
ExactSolution () : Function<dim>(dim+1) {}
virtual void vector_value (const Point<dim> &p,
Vector<double> &value) const;
virtual void vector_gradient (const Point<dim> &p,
std::vector<Tensor<1,dim> > &gradients) const;
};
template <int dim>
void
ExactSolution<dim>::vector_value (const Point<dim> &p,
Vector<double> &values) const
{
Assert (values.size() == dim+1,
ExcDimensionMismatch (values.size(), dim+1));
if (p[dim-1]>=0.5){
values(0) = 200 * pow(p[0],2) * pow(1 - p[0], 2) * p[1] * (1 - p[1]) * (1 - 2*p[1]);
values(1) = -200 * pow(p[1],2) * pow(1 - p[1], 2) * p[0] * (1 - p[0]) * (1 - 2*p[0]);
//values(2) = 10 * pow(p[0] - 0.5, 3) * pow(p[1], 2) + pow(1 - p[0], 3) * pow(p[1] - 0.5, 3) ;
//values(2) = p[0] * (p[0] - 0.5) * (p[0] - 1.0) * (pow(p[1],2) -p[1] + (double) (((double )1.0)/6.0));
//values(2) = p[0] * (p[0] - 0.5) * (p[0] - 1.0) * (pow(p[1],2) -p[1] + (double) (((double )1.0)/6.0));
values(2) = 0;// p[0] * (p[0] - 0.5) * (p[0] - 1.0) * (pow(p[1],3) -1.5*pow(p[1],2) +0.625*p[1] + 0.0625); // cubic L20 pressure
}
else {
values(0) = 200 * pow(p[0],2) * pow(1 - p[0], 2) * p[1] * (1 - p[1]) * (1 - 2*p[1]);
values(1) = -200 * pow(p[1],2) * pow(1 - p[1], 2) * p[0] * (1 - p[0]) * (1 - 2*p[0]);
values(2) = 0.0;
}
}
///////////
//////////Stuff for the schur solver
//// Inverse matrix stuff like that
template <class Matrix, class Preconditioner>
class InverseMatrix : public Subscriptor
{
public:
InverseMatrix (const Matrix &m,
const Preconditioner &preconditioner);
void vmult (Vector<double> &dst,
const Vector<double> &src) const;
private:
const SmartPointer<const Matrix> matrix;
const SmartPointer<const Preconditioner> preconditioner;
};
template <class Matrix, class Preconditioner>
InverseMatrix<Matrix,Preconditioner>::InverseMatrix (const Matrix &m,
const Preconditioner &preconditioner)
:
matrix (&m),
preconditioner (&preconditioner)
{}
template <class Matrix, class Preconditioner>
void InverseMatrix<Matrix,Preconditioner>::vmult (Vector<double> &dst,
const Vector<double> &src) const
{
SolverControl solver_control (src.size(), 1e-6*src.l2_norm());
SolverCG<> cg (solver_control);
dst = 0;
cg.solve (*matrix, dst, src, *preconditioner);
}
template <class Preconditioner>
class SchurComplement : public Subscriptor
{
public:
SchurComplement (const BlockSparseMatrix<double> &system_matrix,
const InverseMatrix<SparseMatrix<double>, Preconditioner> &A_inverse);
void vmult (Vector<double> &dst,
const Vector<double> &src) const;
private:
const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
const SmartPointer<const InverseMatrix<SparseMatrix<double>, Preconditioner> > A_inverse;
mutable Vector<double> tmp1, tmp2;
};
template <class Preconditioner>
SchurComplement<Preconditioner>::
SchurComplement (const BlockSparseMatrix<double> &system_matrix,
const InverseMatrix<SparseMatrix<double>,Preconditioner> &A_inverse)
:
system_matrix (&system_matrix),
A_inverse (&A_inverse),
tmp1 (system_matrix.block(0,0).m()),
tmp2 (system_matrix.block(0,0).m())
{}
template <class Preconditioner>
void SchurComplement<Preconditioner>::vmult (Vector<double> &dst,
const Vector<double> &src) const
{
system_matrix->block(0,1).vmult (tmp1, src);
A_inverse->vmult (tmp2, tmp1);
system_matrix->block(1,0).vmult (dst, tmp2);
}
/////// End of stuff for schur solver
/////////
template <int dim>
void
ExactSolution<dim>::vector_gradient (const Point<dim> &p,
std::vector<Tensor<1,dim> > &gradients) const
{
Assert (gradients.size() == dim+1,
ExcDimensionMismatch (gradients.size(), dim+1));
//
// values(0) = 200 * pow(p[0],2) * pow(1 - p[0], 2) * p[1] * (1 - p[1]) * (1 - 2*p[1]);
// values(1) = -200 * pow(p[1],2) * pow(1 - p[1], 2) * p[0] * (1 - p[0]) * (1 - 2*p[0]);
// values(2) = 10 * pow(p[0] - 0.5, 3) * pow(p[1], 2) + pow(1 - p[0], 3) * pow(p[1] - 0.5, 2);
Tensor<1,dim> g;
g[0] = 200 * ( 2 * p[0] - 6*pow(p[0], 2) + 4 * pow(p[0],3) ) * p[1] * (1 - p[1]) * (1 - 2*p[1]) ;
g[1] = 200 * pow(p[0], 2) * pow(1 - p[0], 2) * (1 - 6*p[1] + 6 * pow(p[1], 2));
gradients[0] = g;
g[0] = -200 * pow(p[1],2) * pow(1 - p[1], 2) * (1 - 6*p[0] + 6 * pow(p[0], 2) );
g[1] = -200 * ( 2 * p[1] - 6*pow(p[1], 2) + 4 * pow(p[1],3) ) * p[0] * (1 - p[0]) * (1 - 2 * p[0]) ;
gradients[1] = g;
//
g[0] = 0;
g[1] = 0;
gradients[2] = g;
}
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+2), dim,
FE_Q<dim>(stokes_degree+1), 1),
laplace_fe (FE_Q<dim>(stokes_degree+2), 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.0)
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());
}
}
// version of setp 46 without schur complement
// 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());
// }
// version with schur complement
template <int dim>
void
FluidStructureProblem<dim>::setup_dofs ()
{
A_preconditioner.reset ();
system_matrix.clear ();
set_active_fe_indices ();
dof_handler.distribute_dofs (fe_collection);
DoFRenumbering::Cuthill_McKee (dof_handler);
std::vector<unsigned int> block_component (dim+1,0);
block_component[dim] = 1;
DoFRenumbering::component_wise (dof_handler, block_component);
{
constraints.clear ();
const FEValuesExtractors::Vector velocities(0);
DoFTools::make_hanging_node_constraints (dof_handler,
constraints);
VectorTools::interpolate_boundary_values (dof_handler,
1,
StokesBoundaryValues<dim>(),
constraints,
fe_collection.component_mask(velocities));
}
// {
// 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::vector<types::global_dof_index> dofs_per_block (2);
DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component);
const unsigned int n_u = dofs_per_block[0],
n_p = dofs_per_block[1];
std::cout << " Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl
<< " Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< " (" << n_u << '+' << n_p << ')'
<< std::endl;
{
BlockDynamicSparsityPattern dsp (2,2);
dsp.block(0,0).reinit (n_u, n_u);
dsp.block(1,0).reinit (n_p, n_u);
dsp.block(0,1).reinit (n_u, n_p);
dsp.block(1,1).reinit (n_p, n_p);
dsp.collect_sizes();
DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false);
sparsity_pattern.copy_from (dsp);
// 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());
solution.reinit (2);
solution.block(0).reinit (n_u);
solution.block(1).reinit (n_p);
solution.collect_sizes ();
system_rhs.reinit (2);
system_rhs.block(0).reinit (n_u);
system_rhs.block(1).reinit (n_p);
system_rhs.collect_sizes ();
}
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]
- stokes_phi_p[i] * stokes_phi_p[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);
}
}
}
std::cout << " Computing preconditioner..." << std::endl << std::flush;
A_preconditioner
= std_cxx11::shared_ptr<typename InnerPreconditioner<dim>::type>(new typename InnerPreconditioner<dim>::type());
A_preconditioner->initialize (system_matrix.block(0,0),
typename InnerPreconditioner<dim>::type::AdditionalData());
}
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;
// -((
// (2*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);
const InverseMatrix<SparseMatrix<double>,
typename InnerPreconditioner<dim>::type>
A_inverse (system_matrix.block(0,0), *A_preconditioner);
Vector<double> tmp (solution.block(0).size());
{
Vector<double> schur_rhs (solution.block(1).size());
A_inverse.vmult (tmp, system_rhs.block(0));
system_matrix.block(1,0).vmult (schur_rhs, tmp);
schur_rhs -= system_rhs.block(1);
SchurComplement<typename InnerPreconditioner<dim>::type>
schur_complement (system_matrix, A_inverse);
SolverControl solver_control (solution.block(1).size(),
1e-6*schur_rhs.l2_norm());
SolverCG<> cg (solver_control);
SparseILU<double> preconditioner;
preconditioner.initialize (system_matrix.block(1,1),
SparseILU<double>::AdditionalData());
InverseMatrix<SparseMatrix<double>,SparseILU<double> >
m_inverse (system_matrix.block(1,1), preconditioner);
cg.solve (schur_complement, solution.block(1), schur_rhs,
m_inverse);
constraints.distribute (solution);
std::cout << " "
<< solver_control.last_step()
<< " outer CG Schur complement iterations for pressure"
<< std::endl;
}
{
system_matrix.block(0,1).vmult (tmp, solution.block(1));
tmp *= -1;
tmp += system_rhs.block(0);
A_inverse.vmult (solution.block(0), tmp);
constraints.distribute (solution);
}
}
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);
ExactSolution<dim> exact_solution;
// QTrapez<1> q_trapez;
// QIterated<dim> quadrature (q_trapez, stokes_degree+2);
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);
Vector<double> cellwise_errors (triangulation.n_active_cells());
Vector<double> lines (solution.size());
VectorTools::integrate_difference (dof_handler, solution, exact_solution,
cellwise_errors, q_collection,
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, q_collection,
// 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, q_collection,
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, q_collection,
VectorTools::Hdiv_seminorm,
&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 << "Errors: ||e_p||_L2 = " << p_l2_error
<< ", ||u_h||_Hdiv = " << u_hdiv
// << " type of dof_handler "
// << typeid(dof_handler).name()
// << " type of solution "
// << typeid(solution).name()
// << " size of solution "
// << solution.size()
// << " type of cellwise_errors "
// << typeid(cellwise_errors).name()
// << " type of quadrature "
// << typeid(quadrature).name()
// << " type of &velocity_mask "
// << typeid(&velocity_mask).name()
<< 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||_{L^2}$", u_l2_error);
convergence_table.add_value("$||u-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);
std::cout << std::scientific<< std::endl;
convergence_table.write_text(std::cout<< std::scientific);
std::string error_filename = "/Users/gs1511/dealii_out/schur_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/schur_gs_error-adaptive.text";
std::ofstream error_table_file2(error_filename2.c_str());
convergence_table.write_text(error_table_file2);
}
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/schur_combi-"
<< 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
*= 1. / stokes_estimated_error_per_cell.l2_norm();
laplace_estimated_error_per_cell
*= 1. / 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 (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)->center()[dim-1] == 0)
cell->face(f)->set_all_boundary_ids(1);
// We then apply an initial refinement before solving for the first
// time. In 3D, there are going to be more degrees of freedom, so we
// refine less there:
triangulation.refine_global (4-dim);
for (unsigned int refinement_cycle = 0; refinement_cycle<10-2*dim-1;
++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;
}
}
}
//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;
//}
////
//
//
//
//
