On 08/19/2016 01:15 PM, thomas stephens wrote:
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.
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.
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.
Well, is system_rhs nonzero?
Best
W.
--
------------------------------------------------------------------------
Wolfgang Bangerth email: [email protected]
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.
