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.
On Thursday, August 18, 2016 at 12:24:23 PM UTC-4, Wolfgang Bangerth wrote:
>
> On 08/18/2016 09:39 AM, thomas stephens wrote:
> > id_X looks like this:
> >
> > |
> > template<intspacedim>
> > classIdentity:publicFunction<spacedim>
> > {
> > 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,
> > 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
> 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).
>
> Best
> W.
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: [email protected]
> <javascript:>
> www: http://www.math.colostate.edu/~bangerth/
>
--
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.
