weakly enforcing an identity - vector mean curvature as the surface laplacian of the identity on a codim-1 manifold
158 views
Skip to first unread message
thomas stephens
Aug 17, 2016, 12:22:07 PM8/17/16
to deal.II User Group
I am working toward simulating Helfrich flow (squared mean curvature flow) on a parametric surface based on [1] Barrett Garcke, and Nurnberg's 2008 SIAM paper Parametric Approximation of Willmore Flow and Related Geoemtric Evolution Equations (pdf attached).
My strategy is to parameterize a Manifold X, and then compute derivatives of the identity on X in order to arrive at normal vectors and curvature information. I know that this road is fraught with difficulty - in particular the normal vector field on a triangulation is not defined at the vertices, and other unpleasantries abound, but for now I plan to see how far I can get.
On to the question: It is a fact that the vector mean curvature k_bar is equal to the surface laplacian of the identity function id_X on the manifold X, k_bar = laplace_beltrami id_X. Equation (4.4) in [1] is the weak form of this identity: (k_bar, eta_bar) = - (surface gradient of id_X, surface gradient of eta_bar) for all eta_bar in some vector-valued finite element space. (For clarity, k_bar = k*nu, where nu is a unit normal vector, and k = k_1 + k_2 is the 'mean' curvature).
I am trying to use dealii to compute k_bar from this expression, but it is not clear to me how set up this problem. It does not look like a PDE to me since the unknown is k_bar, not id_X. If I were to write the standard Poisson Equation: laplace u = f, then in this case, my unknown is f, not u! Maybe I have been staring at this for too long - is there another way to write equation (4.4) to make it appear more familiar to me? Should I really try to solve u = inverse laplace f? (where, in my case, that would be id_X = inverse laplace k_bar). If this is the case, how does one pull off the inverse laplacian in dealii?
What I have so far is an Ellipsoid<dim,spacedim> class (derived from SphericalManifold - thank you Luca Hentai), and a vector-valued Function<spacedim> which is just the identity. My fe_values come from an FESystem<dim,spacedim>, and I'm hoping to be able to use fe_values.shape_grad() or something similar to get at these surface gradients.
Thank you for your time
Wolfgang Bangerth
Aug 17, 2016, 12:48:05 PM8/17/16
to dea...@googlegroups.com
On 08/17/2016 10:22 AM, thomas stephens wrote:
> *On to the question: *It is a fact that the vector mean curvature k_bar
evaluate the curvate at each point in isolation without having to know
anything about the curvate anywhere else.
But it's an equation nonetheless that is best solved in variational
form. Let's say, you approximate k_bar by k_{bar,h} from some finite
element space, then it satisfies the variational equality
(v_h, k_{bar,h})_X = - (nabla_X v_h, nabla_X id_X)_X
On the left hand side, this results in a mass matrix that you will have
to invert to obtain the mean curvature.
Best
W.
--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/
> *On to the question: *It is a fact that the vector mean curvature k_bar
> is equal to the surface laplacian of the identity function id_X on the
> manifold X, k_bar = laplace_beltrami id_X. Equation (4.4) in [1] is
> the weak form of this identity: (k_bar, eta_bar) = - (surface gradient
> of id_X, surface gradient of eta_bar) for all eta_bar in some
> vector-valued finite element space. (For clarity, k_bar = k*nu, where
> nu is a unit normal vector, and k = k_1 + k_2 is the 'mean' curvature).
>
> I am trying to use dealii to compute k_bar from this expression, but it
> is not clear to me how set up this problem. It does not look like a PDE
> to me since the unknown is k_bar, not id_X.
It isn't a PDE because (on an analytically known surface X) you can
> manifold X, k_bar = laplace_beltrami id_X. Equation (4.4) in [1] is
> the weak form of this identity: (k_bar, eta_bar) = - (surface gradient
> of id_X, surface gradient of eta_bar) for all eta_bar in some
> vector-valued finite element space. (For clarity, k_bar = k*nu, where
> nu is a unit normal vector, and k = k_1 + k_2 is the 'mean' curvature).
>
> I am trying to use dealii to compute k_bar from this expression, but it
> is not clear to me how set up this problem. It does not look like a PDE
> to me since the unknown is k_bar, not id_X.
evaluate the curvate at each point in isolation without having to know
anything about the curvate anywhere else.
But it's an equation nonetheless that is best solved in variational
form. Let's say, you approximate k_bar by k_{bar,h} from some finite
element space, then it satisfies the variational equality
(v_h, k_{bar,h})_X = - (nabla_X v_h, nabla_X id_X)_X
On the left hand side, this results in a mass matrix that you will have
to invert to obtain the mean curvature.
Best
W.
--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/
thomas stephens
Aug 18, 2016, 11:39:47 AM8/18/16
to deal.II User Group, bang...@tamu.edu
id_X looks like this:
template<int spacedim>
class Identity: public Function<spacedim>
{
public:
Identity() : Function<spacedim>() {}
virtual void vector_value(const Point<spacedim> &p, Vector<double> &value) const;
virtual double value(const Point<spacedim> &p, const unsigned int component = 0) const;
};
Identity<spacedim> id_X();Do I need to implement nabla_X id_X from W.'s expression (v_h, k_{bar,h})_X = - (nabla_X v_h, nabla_X id_X)_X above, or is there a way to convince dealii to do this for me?
At the moment, the Identity class doesn't know anything about the Manifold object that I have attached to my Triangulation, so at a minimum I figure I would either have to tell Identity about that or attach the Identity Function to something that knows about the manifold, but somehow it doesn't seem right to attach it to another finite element space.
Wolfgang Bangerth
Aug 18, 2016, 12:24:23 PM8/18/16
to thomas stephens, deal.II User Group
On 08/18/2016 09:39 AM, thomas stephens wrote:
> id_X looks like this:
>
> |
> template<intspacedim>
> classIdentity:publicFunction<spacedim>
> id_X looks like this:
>
> |
> template<intspacedim>
> {
> public:
> Identity():Function<spacedim>(){}
>
>
> virtualvoidvector_value(constPoint<spacedim>&p,Vector<double>&value)const;
> virtualdoublevalue(constPoint<spacedim>&p,constunsignedintcomponent
> =0)const;
> };
>
> Identity<spacedim>id_X();
> |
>
>
> Do I need to implement nabla_X id_Xfrom W.'s expression (v_h,
> public:
> Identity():Function<spacedim>(){}
>
>
> virtualvoidvector_value(constPoint<spacedim>&p,Vector<double>&value)const;
> virtualdoublevalue(constPoint<spacedim>&p,constunsignedintcomponent
> =0)const;
> };
>
> Identity<spacedim>id_X();
> |
>
>
> k_{bar,h})_X = - (nabla_X v_h, nabla_X id_X)_X above, or is there a
> way to convince dealii to do this for me?
You can interpolate this function onto a finite element function, and
> way to convince dealii to do this for me?
then compute its gradient.
I suspect that it would be simpler to use the form of the surface
gradient that involves the normal vector. You know what
grad id_X
is (namely the identity matrix), from which you can compute
grad_X id_X
if you know the normal vector (which you should be able to get from X
somehow).
thomas stephens
Aug 19, 2016, 3:15:43 PM8/19/16
to deal.II User Group, tdste...@gmail.com, bang...@colostate.edu
Wolfgang, I appreciate your assistance so far.
I'm getting closer, but I still need some help. I would like to look at this problem as it appears in [1], Equation 27 on p. 11, as the bilinear forms are written as integrals, making things more clear.
Below I have set up what I think to be the system for (k_{bar,h})_X = - (nabla_X v_h, nabla_X id_X)_X - which, in light of Eqn 27 in [1], is really
\integral k_bar_h \dot v_h = \integral nabla_X id_X : nabla_X v_h
typename DoFHandler<dim,spacedim>::active_cell_iterator cell = dof_handler.begin_active(); for (; cell!=dof_handler.end(); ++cell) { cell_mass_matrix = 0; cell_rhs = 0; fe_values.reinit (cell);
// create lhs mass matrix, (k*nu, v)_ij for (unsigned int i=0; i<dofs_per_cell; ++i) for (unsigned int j=0; j<dofs_per_cell; ++j) for (unsigned int q_point=0; q_point<n_q_points; ++q_point) cell_mass_matrix(i,j) += fe_values.shape_value(i,q_point) * fe_values.shape_value(j,q_point) * fe_values.JxW(q_point);
// create rhs vector, (nabla_X id_X, nabla_X v)_i := \int nabla_X id_X : nabla_X v for (unsigned int i=0; i<dofs_per_cell; ++i) for (unsigned int q_point=0; q_point<n_q_points; ++q_point) cell_rhs(i) += (fe_values.shape_grad_component(i,q_point,0)* identity_on_manifold.shape_grad_component(fe_values.normal_vector(q_point),0) + fe_values.shape_grad_component(i,q_point,1)* identity_on_manifold.shape_grad_component(fe_values.normal_vector(q_point),1) + fe_values.shape_grad_component(i,q_point,2)* identity_on_manifold.shape_grad_component(fe_values.normal_vector(q_point),2))* fe_values.JxW(q_point);
cell->get_dof_indices (local_dof_indices); for (unsigned int i=0; i<dofs_per_cell; ++i) { system_rhs(local_dof_indices[i]) += cell_rhs(i);
for (unsigned int j=0; j<dofs_per_cell; ++j) mass_matrix.add (local_dof_indices[i], local_dof_indices[j], cell_mass_matrix(i,j)); } }
Already I am confused by this, since the double contraction : (which is implemented long-hand above) ultimately puts a scalar into cell_rhs(i), whereas I think I am looking for each element of cell_rhs(i) to have spacedim=3 components. I am quite confused about this actually. (The reason why I think elements at cell_rhs(i) should be vector-valued is because each element in k_bar is vector-valued, so in mass_matrix*k_bar = system_rhs, the components of system_rhs should be vectors.
If the above is correct, then I need help getting the mean curvature vector out of this system. If it is the case that I have the pieces to write mass_matrix*mean_curvature_vector = system_rhs, and I need to invert mass_matrix, then how to extract mean_curvature_vector? The following code yields all zeros for mean_curvature_vector:
/// solve mean_curvature_vector = inverse_mass_matrix * system_rhs SolverControl solver_control (mean_curvature_vector.size(), 1e-7 ); SolverCG<> cg_solver (solver_control);
const auto inverse_mass_matrix = inverse_operator(linear_operator(mass_matrix), cg_solver, PreconditionIdentity()); mean_curvature_vector.reinit(dof_handler.n_dofs()); mean_curvature_squared.reinit(dof_handler.n_dofs()); system_rhs *= -1.0; inverse_mass_matrix.vmult(mean_curvature_vector,system_rhs);
std::cout << "mean curvature vector" << mean_curvature_vector << std::endl; mean_curvature_squared = mean_curvature_vector*mean_curvature_vector; std::cout << "mean curvature squared " << mean_curvature_squared << std::endl;
That is, I do not think I am handling these data structures correctly - perhaps that LinearOperator inverse_mass_matrix is not just an ordinary dealii::Matrix that I can then apply to dealii::Vectors.
thomas stephens
Aug 19, 2016, 3:19:26 PM8/19/16
to deal.II User Group, tdste...@gmail.com, bang...@colostate.edu
(continued from above) The reference [1] above is Parametric FEM for Geometric Biomembranes, Nochetto, Bonito, Pauletti 2011 preprint (pdf attached). My current .cc file is attached as well.
Wolfgang Bangerth
Aug 22, 2016, 10:12:48 AM8/22/16
to deal.II user group
leads to just a scalar equation on both the left and right hand side. Take a
look how this works for other vector-valued problems in the documentation
module on vector-valued problems.
In your code above, it is not clear to me what 'identity_in_manifold' is. Does
this simply represent X? If so, aren't you just computing
grad phi_i * grad X
? You need to somehow use the grad_Gamma instead, though.
> If the above is correct, then I need help getting the mean curvature vector
> out of this system. If it is the case that I have the pieces to write
> mass_matrix*mean_curvature_vector = system_rhs, and I need to invert
> mass_matrix, then how to extract mean_curvature_vector? The following code
> yields all zeros for mean_curvature_vector:
>
> |
> /// solve mean_curvature_vector = inverse_mass_matrix * system_rhs
> SolverControl solver_control (mean_curvature_vector.size(), 1e-7 );
> SolverCG<> cg_solver (solver_control);
>
> const auto inverse_mass_matrix = inverse_operator(linear_operator(mass_matrix),
> cg_solver,
> PreconditionIdentity());
> mean_curvature_vector.reinit(dof_handler.n_dofs());
> mean_curvature_squared.reinit(dof_handler.n_dofs());
> system_rhs *= -1.0;
> inverse_mass_matrix.vmult(mean_curvature_vector,system_rhs);
>
>
> std::cout << "mean curvature vector" << mean_curvature_vector << std::endl;
> mean_curvature_squared = mean_curvature_vector*mean_curvature_vector;
> std::cout << "mean curvature squared " << mean_curvature_squared << std::endl;
>
> |
>
>
> That is, I do not think I am handling these data structures correctly -
> perhaps that LinearOperator inverse_mass_matrix is not just an ordinary
> dealii::Matrix that I can then apply to dealii::Vectors.
thomas stephens
Aug 23, 2016, 11:01:03 AM8/23/16
to deal.II User Group, bang...@colostate.edu
I was able to solve for the vector mean curvature using the weak form of the identity k_bar = laplacian_X id_X, where X is a codimension 1 manifold without boundary. The image above is a plot of the square of the mean curvature on an ellipsoid with semi principle axes a=1,b=2,c=3.
For concreteness, X is a 2-d surface embedded in R^3, id_X : (x,y,z) |----> (x,y,z) is simply the identity on X, k_bar is the vector mean curvature, and has the form k*nu, where k = k_1 + k_2 is the scalar curvature, and nu is a normal vector on X. The variational form of the identity is (v_bar, k_bar)_X = - (nabla_X v_bar, nabla_X id_X)_X, for v_bar coming from a vector-valued finite element space which appears below as FE_Q<dim,spacedim> = FE_Q<2,3>. A good starting point for this problem is [1], and the weak identity is printed as Eqn (4.4).
Now the solution I seek is the vector-valued mean curvature k_bar. Following W.'s suggestion above, I built a mass_matrix for the lhs and filled out the rhs with known information as system_rhs, and I want to solve mass_matrix*k_bar = system_rhs for the vector-valued k_bar. I have naively implemented this componentwise below, so the matrix equation is solved three times: solver.solve(mass_matrix, k_bar_x, system_rhs_x), solver.solve(mass_matrix, k_bar_y, system_rhs_y), solver.solve(mass_matrix, k_bar_z, system_rhs_z), and the data I am looking for is the mean curvature squared, so at the bottom of the solve() function I compute the scalar product mean_curvature_squared = k_bar*k_bar.
References [1], [2], and [3] below are very good. I adapted Tutorial step-38 to fit my needs. To get the geometry I trivially extended SphericalManifold to get an Ellipsoid, this .cc file is attached. The identity id_X is just a Function, and appears in the code below.
There are other ways to compute the (vector and scalar) mean curvature, and for my next trick I will attempt to compute the trace of the Weingarten map, nabla_X nu. According to Eqn (4.3) in [1], this will yield the scalar mean curvature k = k_1 + k_2. (Please don't fret - in the literature 'mean curvature' sometimes means the arithmetic mean of the principle curvatures and sometimes means the 'total curvature', just the sum of principle curvatures without the normalization! I am perpetuating the ambiguity by calling k_1 + k_2 the 'mean curvature'.)
Now for the code:
/* ---------------------------------------------------------------------
* mean_curvature_via_step38.cc
*
* MODIFIED VERSION OF STEP-38, TOM STEPHENS, August 2016
*
* Copyright (C) 2010 - 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.
*
* ---------------------------------------------------------------------
*
* Authors: Andrea Bonito, Sebastian Pauletti.
*/
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <stdio.h>
#include <fstream>
#include <iostream>
// additional cpp code that is reusable
#include "../../utilities/my_manifolds_lib.cc"
namespace Step38
{
using namespace dealii;
template <int spacedim>
class LaplaceBeltramiProblem
{
public:
LaplaceBeltramiProblem (const unsigned degree = 2);
void run ();
private:
static const unsigned int dim = spacedim-1;
void make_grid_and_dofs ();
void assemble_system ();
void solve ();
void output_results () const;
void compute_error () const;
Triangulation<dim,spacedim> triangulation;
FE_Q<dim,spacedim> fe;
double degree = 2;
DoFHandler<dim,spacedim> dof_handler;
MappingQ<dim, spacedim> mapping;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution_x;
Vector<double> solution_y;
Vector<double> solution_z;
Vector<double> system_rhs_x;
Vector<double> system_rhs_y;
Vector<double> system_rhs_z;
Vector<double> mean_curvature_squared;
};
template <int spacedim>
class Identity : public Function<spacedim>
{
/*{{{*/
public:
Identity() : Function<spacedim>() {}
virtual void vector_value_list (const std::vector<Point<spacedim> > &points,
std::vector<Vector<double> > &value_list) const;
virtual void vector_value (const Point<spacedim> &p, Vector<double> &value) const;
virtual double value (const Point<spacedim> &p, const unsigned int component = 0) const;
virtual Tensor<2,spacedim> shape_grad(const Tensor<1,spacedim> &unit_normal) const;
virtual Tensor<1,spacedim> shape_grad_component(const Tensor<1,spacedim> &unit_normal, const unsigned int component) const;
/*}}}*/
};
template<int spacedim>
double Identity<spacedim>::value(const Point<spacedim> &p, const unsigned int component) const
{
/*{{{*/
return p(component);
/*}}}*/
}
template<int spacedim>
void Identity<spacedim>::vector_value(const Point<spacedim> &p, Vector<double> &value) const
{
/*{{{*/
for (unsigned int c=0; c<this->n_components; ++c)
{
value(c) = Identity<spacedim>::value(p,c);
}
/*}}}*/
}
template <int spacedim>
void Identity<spacedim>::vector_value_list (const std::vector<Point<spacedim> > &points,
std::vector<Vector<double> > &value_list) const
{
/*{{{*/
Assert (value_list.size() == points.size(),
ExcDimensionMismatch (value_list.size(), points.size()));
const unsigned int n_points = points.size();
for (unsigned int p=0; p<n_points; ++p)
Identity<spacedim>::vector_value (points[p],
value_list[p]);
/*}}}*/
}
template <int spacedim>
Tensor<2,spacedim> Identity<spacedim>::shape_grad(const Tensor<1,spacedim> &unit_normal) const
{
/*{{{*/
Tensor<2,spacedim> eye;
eye = 0; eye[0][0] = 1; eye[1][1] = 1; eye[2][2] = 1;
Tensor<2,spacedim> nnT;
nnT = outer_product(unit_normal,unit_normal);
return eye - nnT;
/*}}}*/
}
template <int spacedim>
Tensor<1,spacedim> Identity<spacedim>::shape_grad_component(const Tensor<1,spacedim> &unit_normal, const unsigned int component) const
{
/*{{{*/
Tensor<2,spacedim> full_shape_grad = shape_grad(unit_normal);
Tensor<1,spacedim> grad_component;
grad_component[0] = full_shape_grad[component][0];
grad_component[1] = full_shape_grad[component][1];
grad_component[2] = full_shape_grad[component][2];
return grad_component;
/*}}}*/
}
template <int spacedim>
LaplaceBeltramiProblem<spacedim>::LaplaceBeltramiProblem (const unsigned degree)
:
fe(degree),
dof_handler(triangulation),
mapping (degree)
{}
template <int spacedim>
void LaplaceBeltramiProblem<spacedim>::make_grid_and_dofs ()
{
double a = 1; double b = 2; double c = 3;
Point<spacedim> center(0,0,0);
static Ellipsoid<dim,spacedim> ellipsoid(a,b,c,center);
GridGenerator::hyper_sphere(triangulation,center, 1);
triangulation.set_all_manifold_ids(0);
GridTools::transform(std_cxx11::bind(&Ellipsoid<dim,spacedim>::grid_transform, &ellipsoid, std_cxx11::_1),
triangulation);
triangulation.set_manifold (0, ellipsoid);
triangulation.refine_global(4);
std::cout << "Surface mesh has " << triangulation.n_active_cells()
<< " cells."
<< std::endl;
dof_handler.distribute_dofs (fe);
std::cout << "Surface mesh has " << dof_handler.n_dofs()
<< " degrees of freedom."
<< std::endl;
DynamicSparsityPattern dsp (dof_handler.n_dofs(), dof_handler.n_dofs());
DoFTools::make_sparsity_pattern (dof_handler, dsp);
sparsity_pattern.copy_from (dsp);
system_matrix.reinit (sparsity_pattern);
solution_x.reinit (dof_handler.n_dofs());
solution_y.reinit (dof_handler.n_dofs());
solution_z.reinit (dof_handler.n_dofs());
system_rhs_x.reinit (dof_handler.n_dofs());
system_rhs_y.reinit (dof_handler.n_dofs());
system_rhs_z.reinit (dof_handler.n_dofs());
}
template <int spacedim>
void LaplaceBeltramiProblem<spacedim>::assemble_system ()
{
system_matrix = 0;
system_rhs_x = 0;
system_rhs_y = 0;
system_rhs_z = 0;
Identity<spacedim> identity_on_manifold;
const QGauss<dim> quadrature_formula(2*fe.degree);
FEValues<dim,spacedim> fe_values (mapping, fe, quadrature_formula,
update_values |
update_normal_vectors |
update_gradients |
update_quadrature_points |
update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs_x (dofs_per_cell);
Vector<double> cell_rhs_y (dofs_per_cell);
Vector<double> cell_rhs_z (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
for (typename DoFHandler<dim,spacedim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
cell!=endc; ++cell)
{
cell_matrix = 0;
cell_rhs_x = 0;
cell_rhs_y = 0;
cell_rhs_z = 0;
fe_values.reinit (cell);
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
{
cell_matrix(i,j) += fe_values.shape_value(i,q_point) *
fe_values.shape_value(j,q_point) *
fe_values.JxW(q_point);
}
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
{
cell_rhs_x(i) += fe_values.shape_grad(i,q_point)*
identity_on_manifold.shape_grad_component(fe_values.normal_vector(q_point),0)*
fe_values.JxW(q_point);
cell_rhs_y(i) += fe_values.shape_grad(i,q_point)*
identity_on_manifold.shape_grad_component(fe_values.normal_vector(q_point),1)*
fe_values.JxW(q_point);
cell_rhs_z(i) += fe_values.shape_grad(i,q_point)*
identity_on_manifold.shape_grad_component(fe_values.normal_vector(q_point),2)*
fe_values.JxW(q_point);
}
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
system_matrix.add (local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i,j));
system_rhs_x(local_dof_indices[i]) += cell_rhs_x(i);
system_rhs_y(local_dof_indices[i]) += cell_rhs_y(i);
system_rhs_z(local_dof_indices[i]) += cell_rhs_z(i);
}
}
}
template <int spacedim>
void LaplaceBeltramiProblem<spacedim>::solve ()
{
SolverControl solver_control (solution_x.size(), 1e-12 );
SolverCG<> cg (solver_control);
SolverGMRES<> gmres (solver_control);
//cg.solve (system_matrix, solution_x, system_rhs_x, preconditioner);
cg.solve (system_matrix, solution_x, system_rhs_x, PreconditionIdentity());
std::cout << "Solved x component" << std::endl;
cg.solve (system_matrix, solution_y, system_rhs_y, PreconditionIdentity());
std::cout << "Solved y component" << std::endl;
cg.solve (system_matrix, solution_z, system_rhs_z, PreconditionIdentity());
std::cout << "Solved z component" << std::endl;
mean_curvature_squared.reinit (dof_handler.n_dofs());
double avg = 0; double summ = 0;
for (unsigned int i=0; i<dof_handler.n_dofs(); ++i )
{
mean_curvature_squared(i) = pow(solution_x(i),2) + pow(solution_y(i),2) + pow(solution_z(i),2);
summ += mean_curvature_squared(i);
}
avg = summ/dof_handler.n_dofs();
std::cout << "avg mean curvature: " << avg << std::endl;
}
template <int spacedim>
void LaplaceBeltramiProblem<spacedim>::output_results () const
{
DataOut<dim,DoFHandler<dim,spacedim> > data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (mean_curvature_squared,
"solution",
DataOut<dim,DoFHandler<dim,spacedim> >::type_dof_data);
data_out.build_patches (mapping,
mapping.get_degree());
std::string filename ("solution-");
filename += static_cast<char>('0'+spacedim);
filename += "d.vtk";
std::ofstream output (filename.c_str());
data_out.write_vtk (output);
}
template <int spacedim>
void LaplaceBeltramiProblem<spacedim>::run ()
{
make_grid_and_dofs();
std::cout << "grid and dofs made " << std::endl;
assemble_system ();
std::cout << "system assembled " << std::endl;
solve ();
std::cout << "solved " << std::endl;
output_results ();
std::cout << "results written" << std::endl;
}
}
int main ()
{
try
{
using namespace dealii;
using namespace Step38;
LaplaceBeltramiProblem<3> laplace_beltrami;
laplace_beltrami.run();
return 0;
}
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::cout << "problem! " << std::endl;
return 2;
}
}
Thank you Wolfgang for your help.
> Already I am confused by this, since the double contraction : (which is
> implemented long-hand above) ultimately puts a scalar into cell_rhs(i),
> whereas I think I am looking for each element of cell_rhs(i) to have
> spacedim=3 components. I am quite confused about this actually. (The reason
> why I think elements at cell_rhs(i) should be vector-valued is because each
> element in k_bar is vector-valued, so in mass_matrix*k_bar = system_rhs, the
> components of system_rhs should be vectors.
No, that doesn't make any sense. For each i, the weak form of the equation
leads to just a scalar equation on both the left and right hand side. Take a
look how this works for other vector-valued problems in the documentation
module on vector-valued problems.
You're right, none of that earlier post made any sense. That's not to say that this post makes sense, either, but that last one was pretty bad. I still think the solution should be vector-valued, and I will now try to make this look less silly by using vector-valued elements rather than solving this componentwise.
In your code above, it is not clear to me what 'identity_in_manifold' is. Does
this simply represent X? If so, aren't you just computing
grad phi_i * grad X
? You need to somehow use the grad_Gamma instead, though.
I think my Identity function implements the shape_grad() correctly here. I used the normal vector from the finite elements: fe_values.normal_vector(q_point) because I ultimately want to evolve this surface and I am hoping that the fe_values will continue to hold the normals I am looking for. This idea may come crashing down.
> If the above is correct, then I need help getting the mean curvature vector
> out of this system. If it is the case that I have the pieces to write
> mass_matrix*mean_curvature_vector = system_rhs, and I need to invert
> mass_matrix, then how to extract mean_curvature_vector? The following code
> yields all zeros for mean_curvature_vector:
>
> |
> /// solve mean_curvature_vector = inverse_mass_matrix * system_rhs
> SolverControl solver_control (mean_curvature_vector.size(), 1e-7 );
> SolverCG<> cg_solver (solver_control);
>
> const auto inverse_mass_matrix = inverse_operator(linear_operator(mass_matrix),
> cg_solver,
> PreconditionIdentity());
> mean_curvature_vector.reinit(dof_handler.n_dofs());
> mean_curvature_squared.reinit(dof_handler.n_dofs());
> system_rhs *= -1.0;
> inverse_mass_matrix.vmult(mean_curvature_vector,system_rhs);
>
>
> std::cout << "mean curvature vector" << mean_curvature_vector << std::endl;
> mean_curvature_squared = mean_curvature_vector*mean_curvature_vector;
> std::cout << "mean curvature squared " << mean_curvature_squared << std::endl;
>
> |
>
>
> That is, I do not think I am handling these data structures correctly -
> perhaps that LinearOperator inverse_mass_matrix is not just an ordinary
> dealii::Matrix that I can then apply to dealii::Vectors.
Well, is system_rhs nonzero?
It doesn't matter what system_rhs was in that previous implementation - the whole thing was 'garbage in, garbage out'. Thanks for being a patient sounding board.
[1] Barrett, Garcke, Nurnberg, Parametric Approximation of Willmore Flow and Related Geometric Evolution Equations, SIAM 2008 (pdf attached at top of this thread)
[2] Soren Bartels, Numerical Methods for Nonlinear PDE, Springer 2015 (this is a book, and contains matlab code to solve this and similar problems)
[3] Bonito, Nochetto, Pauletti, Parametric FEM for Geometric Biomembranes, 2011 preprint (pdf attached, also see past and present work by the authors, and see references for Dziuk and Elliott's work)
Wolfgang Bangerth
Aug 23, 2016, 12:51:53 PM8/23/16
to thomas stephens, deal.II User Group
Thomas,
> I was able to solve for the vector mean curvature using the weak form of
> the identity k_bar = laplacian_X id_X, where X is a codimension 1
> manifold without boundary. The image above is a plot of the square of
> the mean curvature on an ellipsoid with semi principle axes a=1,b=2,c=3.
numerical values are correct as well?
I think this would actually make for a really nice addition to the set
of tutorials, or the code gallery at
http://dealii.org/code-gallery.html
Would you be interested in getting it there? (Maybe after rewriting it
in such a way that you solve for all three components at once?)
Cheers
thomas stephens
Aug 23, 2016, 5:44:23 PM8/23/16
to deal.II User Group, tdste...@gmail.com, bang...@colostate.edu
On Tuesday, August 23, 2016 at 12:51:53 PM UTC-4, Wolfgang Bangerth wrote:
Thomas,
> I was able to solve for the vector mean curvature using the weak form of
> the identity k_bar = laplacian_X id_X, where X is a codimension 1
> manifold without boundary. The image above is a plot of the square of
> the mean curvature on an ellipsoid with semi principle axes a=1,b=2,c=3.
Great, this looks roughly correct I would say. Have you checked that the
numerical values are correct as well?
I implemented an exact solution for the mean curvature for ellipsoids (it's a slightly ugly expression, I just took it from here: http://math.stackexchange.com/a/540820), and I computed the L2 and Linfty errors between my computed values and that expression. The result: Convergence wrt mesh refinement - about an order of magnitude decrease in each norm per increment of refine_global(). That's not a tidy way to put it, but it was quick and easy.
I think this would actually make for a really nice addition to the set
of tutorials, or the code gallery at
http://dealii.org/code-gallery.html
Would you be interested in getting it there? (Maybe after rewriting it
in such a way that you solve for all three components at once?)
Yes, I am very interested in doing this ( but it will have to wait until September), thank you for asking.
Wolfgang Bangerth
Aug 23, 2016, 5:55:34 PM8/23/16
to thomas stephens, deal.II User Group
> I implemented an exact solution for the mean curvature for ellipsoids
> (it's a slightly ugly expression, I just took it from
> here: http://math.stackexchange.com/a/540820), and I computed the L2 and
> Linfty errors between my computed values and that expression. The
> result: Convergence wrt mesh refinement - about an order of magnitude
> decrease in each norm per increment of refine_global(). That's not a
> tidy way to put it, but it was quick and easy.
> I think this would actually make for a really nice addition to the set
> of tutorials, or the code gallery at
> http://dealii.org/code-gallery.html
> <http://dealii.org/code-gallery.html>
> Would you be interested in getting it there? (Maybe after rewriting it
> in such a way that you solve for all three components at once?)
>
>
> Yes, I am very interested in doing this ( but it will have to wait until
> September), thank you for asking.
to go forward with this!
Reply all
Reply to author
Forward
0 new messages