Pankaj,
Hello,
> I am trying to write a code to calculate buoyancy force which is the
> product of temperature and gravity. In my case temperature is available.
>
>
>
> TrilinosWrappers::MPI::Vector temperature ;//given
> const dealii::Point<dim> gravity = -( (dim == 2) ? (dealii::Point<dim> (
> 0,1)) :
> (dealii::Point<dim> (0,0,1)) );
>
> TrilinosWrappers::MPI::Vector& Buoyancy_force = constant * temperature *
> gravity ; //This is what i want to do
>
If you want to calculate your buyoancy as TrilinosWrappers::MPI::Vector,
you need a corresponding DoFHandler.
>
>
>
> //below is the idea that i am not sure of
>
> typename DoFHandler<dim>::active_cell_iterator
> cell = fe_velocity_space.dh_begin_active(),
> endc = fe_velocity_space.dh_end(),
> cell_t = fe_temp_space.dh_begin_active();
> for (; cell != endc; ++cell, ++cell_t) {
> if (cell->is_locally_owned() == false)
> continue;
>
> b = 0;
>
> fe_velocity_values.reinit(cell);
> fe_temp_values.reinit(cell_t);
> const FEValuesViews::Vector<dim>& fe_vector_values =
> fe_velocity_values[FEValuesExtractors::Vector(0)];
> fe_temp_values.get_function_values(temperature, temp_cell);
>
> double Beta = 0.04;
>
> for (unsigned int q=0; q<nqcell; ++q){
> const double& temp_q = temp_cell[q];
> const Tensor<1, dim> force = - Beta *(temp_q)*fe_values.
> JxW(q)*gravity ;
> //Need help here
>
>
>
> cell->get_dof_indices(dof_indices);
>
> fe_velocity_space.get_constraints().distribute_local_to_global
> (b, dof_indices, buoyancy_force);
>
> }
>
I guess that in all places, in which you are using fe_velocity, you mean to
use a FiniteElement for the buoyancy force? The variable b should be force,
correct?
What you are doingso far is to compute the buoyancy force in each
quadrature point. If your FiniteElement is interpolating in its
support_points (as FE_Q), then
you want to use as quadrature points in the FEValues objects
fe_velocity_values and fe_temp_values the unit_support_points of the
FiniteElement [1].
The number of quadrature points you get after reinitializing the FEValues
objects equals the number of DoFs in each component multiplied by the
number of components and contains the mapped support points. Then you can
simply assign the values you computed above. In particular, you can do
something like:
Quadrature q(fe.get_unit_support_points());
FEValues fe_values (..., q, update_q_points);
...
for (unsigned int q=0; q<nqcell; ++q){
const double& temp_q = temp_cell[q];
const Tensor<1, dim> force = - Beta *(temp_q)*fe_values.JxW(q)*gravity ;
const unsigned int component =
FiniteElement::system_to_component_index(q).first;
b(q)=force(component)
}
...
Note that in this snippet, the temperature in the same quadrature point is
calculated multiply.
Best,
Daniel
[1]
https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-to-get-the-mapped-position-of-support-points-of-my-element
--
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.
