Dear all,
I am hoping someone can shed a bit more light on a problem I am having. The
background of the system of equations i am solving: I have 4 DAEs (for pf,
pr, vr, vf) and 2 time dependent equations (phi, t) as a system of 6
equations. I am solving the 4daes first then using the solutions into the
time dependent equations (which need the solutions in the DAEs), and I am
using these two solutions to rerun the 4 DAEs and so on.
The right hand side of the weak forms of the DAEs all require information
about each other, so at the initial timestep I am providing a guess, but
from then on, i am using get_function_values or
get_function_gradients/divergences. I have been using these fine for most
of my variables, but i am having an issue with one of them.
My 4 DAEs look basically like the Stokes' system, with 1 DAEs in mixed
formulation (very similar to Darcy's) and another non-differential
algebraic equation for vf. Now, for this equation/variable, I am using:
FESystem<dim> vf_fe;
DoFHandler<dim> vf_dof_handler;
ConstraintMatrix vf_constraints;
SparsityPattern vf_sparsity_pattern;
SparseMatrix<double> vf_system_matrix;
SparseMatrix<double> vf_mass_matrix;
Vector<double> vf_solution;
Vector<double> vf_initial_solution;
Vector<double> vf_system_rhs;
and vf_fe (FE_Q<dim>(pf_degree+1), dim) as vf is vector valued.
For the others, like the example codes (step-20/22/31), i am using
BlockVectors etc. eg,
rock_fe (FE_Q<dim>(pr_degree+1), dim,
FE_Q<dim>(pr_degree), 1)
with BlockSparsityPattern rock_sparsity_pattern;
BlockSparseMatrix<double> rock_system_matrix;
BlockVector<double> rock_solution;
BlockVector<double> rock_system_rhs;
and so on.
I am of course using const FEValuesExtractors::Vector velocities (0); etc
for most of the things I need. now out of the following, in the loop over
all the cells: I have
pf_fe_values[pressure].get_function_values (pf_solution,
pf_values);
pf_fe_values[pressure].get_function_gradients (pf_solution,
grad_pf_values);
phi_fe_values.get_function_values (phi_solution,
phi_values);
phi_fe_values.get_function_gradients (phi_solution,
grad_phi_values);
vf_fe_values.get_function_values (vf_solution, vf_values);
vf_fe_values.get_function_divergences (vf_solution,
div_vf_values);
and the LAST one for get_function_divergences does not seem to work. It is
telling me there is no member function called get_function_divergeneces. I
am a little confused as elsewhere i
use vr_fe_values[velocities].get_function_divergences (rock_solution,
div_vr_values); and this works fine.
I have also tried vf_fe_values[velocities].get_function_divergences
(vf_solution, div_vf_values); - now this does not give me the same error,
the program runs but i get:
An error occurred in line <2364> of file
</scratch/leej/dealnew/dealii-8.5.0/source/fe/fe_values.cc> in function
void dealii::internal::do_function_values(const Number*, const
dealii::Table<2, double>&, const dealii::FiniteElement<dim, spacedim>&,
const std::vector<unsigned int>&,
dealii::VectorSlice<std::vector<VectorType> >&, bool, unsigned int) [with
int dim = 2; int spacedim = 2; VectorType = dealii::Vector<double>; Number
= double]
The violated condition was:
(values[i].size()) == (result_components)
Additional information:
Dimension 0 not equal to 2.
I also get the same no member function call when used with the
get_function_values on the same variable with FEValuesExtractors and adding
[velocities]
I would be grateful if someone could shed a bit more light on why I am
getting these errors. I particularly do not understand why it is giving me
a dimension error either. It seems like quite a trivial thing but it seems
not? I wondered whether it is because it is not aware that the solution is
a vector, but I'm pretty sure the fact that I have an FESystem sorts this
out?
As further information, note that i am using
std::vector<Vector<double> > vf_values (n_q_points);
std::vector<double> div_vf_values (n_q_points);
I have also looped over the right cells etc using
typename DoFHandler<dim>::active_cell_iterator
vf_cell = vf_dof_handler.begin_active();
as well as intialising correctly.
Before this assmebly bit is run, the initial vf_solution is just
interpolated with
VectorTools::interpolate(vf_dof_handler, vfInitialFunction<dim>(),
vf_solution);
where it is just zero
template <int dim>
double
vfInitialFunction<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
if (component == 0)
return 0.0;
else if (component == 1)
return 0.0*p[0];
return 0.0;
}
Many thanks for your help in advance.
--
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.
