I am filling it up like that, and then (when filling system_matrix) I would
like to take a look at the values:
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
{
for(size_t i = 0; i < dofs_per_cell; ++i)
{
gradient_N[i] = fe_values[N_density].gradient(i, q_point
);
gradient_TE[i] = fe_values[E_temperature].gradient(i,
q_point);
gradient_TL[i] = fe_values[L_temperature].gradient(i,
q_point);
value_N[i] = fe_values[N_density].value(i, q_point);
value_TE[i] = fe_values[E_temperature].value(i, q_point
);
value_TL[i] = fe_values[L_temperature].value(i, q_point
);
}
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
std::cout << "N_value = (" << value_N[i] << ",\t " <<
value_N[j] << ")\n, TE_value = (" << value_TE[i] << ",\t " << value_TE[j] <<
")\n,
TL_value = (" << value_TL[i] << ",\t " << value_TL[j] << ")\n";
std::cout << "N_grad = (" << gradient_N[i] << ",\t "
<< gradient_N[j] << ")\n, TE_grad = (" << gradient_TE[i] << ",\t " <<
gradient_TE[j] << ")\n, TL_grad = (" << gradient_TL[i] << ",\t " <<
gradient_TL[j] << ")\n";
getchar();
volatile double source_terms = ( 1/*Alternative to
calculation terms*/*fe_values.JxW(q_point));
cell_matrix(i,j) += source_terms;
volatile double mass_terms = (scalar_product(value_N[i
], value_N[j]) + scalar_product(value_TE[i], value_TE[j])+scalar_product(
value_TL[i], value_TL[j]))*
fe_values.JxW(q_point);
cell_mass_matrix(i,j) +=mass_terms;
}
}
My filling-function is
template <int dim>
class InitialValues : public Function<dim>
{
private:
equation_class eq_class;
physics_equations equations;
public:
InitialValues(const equation_class equ_class) : Function<dim>(3*dim
), eq_class(equ_class) {}
virtual double value (const Point<dim> &p, const unsigned int
component = 0) const;
virtual void vector_value(const Point<dim> &p, Vector<double> &
values) const;
};
template <int dim>
double InitialValues<dim>::value(const Point<dim> &p, const unsigned int
component) const
{
(void) p;
Assert(component < this->n_components, ExcIndexRange(component, 0,
this->n_components));
if(component < dim)
return 1e18;
else
if(component < 2*dim && component >= dim)
return 293;
else
return 293;
}
template <int dim>
void InitialValues<dim>::vector_value(const Point<dim> &p, Vector
<double> &values) const
{
for(size_t i = 0; i < this->n_components; ++i)
values(i) = InitialValues<dim>::value(p, i);
}
and it is called as
VectorTools::project(dof_handler, constraint_matrix, QGauss<dim>(3),
InitialValues<dim>(carrier_density), solution);
Now, when looping over the different values, I get an output of
N_value = (0.324665 0 0, 0.324665 0 0)
, TE_value = (0 0 0, 0 0 0)
, TL_value = (0 0 0, 0 0 0)
N_grad = (-1.20419 -1.20419 -1.20419 0 0 0 0 0 0, -1.20419 -1.20419 -
1.20419 0 0 0 0 0 0)
, TE_grad = (0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0)
, TL_grad = (0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0)
N_value = (0.324665 0 0, 0 0.324665 0)
, TE_value = (0 0 0, 0 0 0)
, TL_value = (0 0 0, 0 0 0)
N_grad = (-1.20419 -1.20419 -1.20419 0 0 0 0 0 0, 0 0 0 -1.20419 -
1.20419 -1.20419 0 0 0)
, TE_grad = (0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0)
, TL_grad = (0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0)
N_value = (0.324665 0 0, 0 0 0.324665)
, TE_value = (0 0 0, 0 0 0)
, TL_value = (0 0 0, 0 0 0)
N_grad = (-1.20419 -1.20419 -1.20419 0 0 0 0 0 0, 0 0 0 0 0 0 -
1.20419 -1.20419 -1.20419)
, TE_grad = (0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0)
, TL_grad = (0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0)
N_value = (0.324665 0 0, 0 0 0)
, TE_value = (0 0 0, 0.324665 0 0)
, TL_value = (0 0 0, 0 0 0)
N_grad = (-1.20419 -1.20419 -1.20419 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0)
, TE_grad = (0 0 0 0 0 0 0 0 0, -1.20419 -1.20419 -1.20419 0 0 0 0 0 0)
, TL_grad = (0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0)
N_value = (0.324665 0 0, 0 0 0)
, TE_value = (0 0 0, 0 0.324665 0)
, TL_value = (0 0 0, 0 0 0)
N_grad = (-1.20419 -1.20419 -1.20419 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0)
, TE_grad = (0 0 0 0 0 0 0 0 0, 0 0 0 -1.20419 -1.20419 -1.20419 0 0 0)
, TL_grad = (0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0)
N_value = (0.324665 0 0, 0 0 0)
, TE_value = (0 0 0, 0 0 0.324665)
, TL_value = (0 0 0, 0 0 0)
N_grad = (-1.20419 -1.20419 -1.20419 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0)
, TE_grad = (0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 -1.20419 -1.20419 -1.20419)
, TL_grad = (0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0)
Why is there only one value not equal to zero, and why is it moving?
Shouldn't all the values be set already?
Am Freitag, 4. August 2017 14:35:07 UTC+2 schrieb Wolfgang Bangerth:
>
> On 08/04/2017 05:41 AM, Maxi Miller wrote:
> > Is an alternative to check the components, and set the conditions
> accordingly?
> > I wrote it like
> > |
> > template<intdim>
> >
> doubleInitialValues<dim>::value(constPoint<dim>&p,constunsignedintcomponent)const
>
>
> > {
> > (void)p;
> > Assert(component
> > <this->n_components,ExcIndexRange(component,0,this->n_components));
> > if(component <dim)
> > returnequations.initial_carrier_density;
> > else
> > if(component <2*dim &&component >=dim)
> > returnequations.ambient_temperature;
> > else
> > returnequations.ambient_temperature;
> > }
> > |
> >
>
> That's how I typically do things.
> 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.
