Hello Professor,
I deleted my question since I thought I had found the answer(in a previous
post of yours) But I guess I am still stuck in the same problem.
What I have done in step 23 is to only make changes in the existing run
function and assemble my lhs matrices and rhs vectors cell wise.
Here is the run function I have tried to implement but the total energy
keeps increasing. Also, since the F values are considered zero, I thought I
should bother myself with F before I get this running. For the past three
days I have been trying to find what I have been missing but, alas :
Thanking you for your time.
template <int dim>
void WaveEquation<dim>::run ()
{
setup_system();
QGauss<dim> quadrature_formula(3);
VectorTools::project (dof_handler, constraints, QGauss<dim>(3),
InitialValuesU<dim>(),
old_solution_u);
VectorTools::project (dof_handler, constraints, QGauss<dim>(3),
InitialValuesV<dim>(),
old_solution_v);
Vector<double> tmp (solution_u.size());
Vector<double> forcing_terms (solution_u.size());
//***********************************************************************************************************************
// Code Substituted
//***********************************************************************************************************************
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
//***********************************************************************************************************************
std::cout << "Dofs per cell : " << dofs_per_cell << std::endl;
std::cout << "Quadrature formula size : " << n_q_points << std::endl;
for (; time<=5; time+=time_step, ++timestep_number)
{
std::cout << "Time step " << timestep_number
<< " at t=" << time
<< std::endl;
//***********************************************************************************************************************
// Code Substituted
//***********************************************************************************************************************
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
Vector<double> cell_rhs_1(dofs_per_cell);
Vector<double> cell_rhs_2(dofs_per_cell);
std::vector<types::global_dof_index>
local_dof_indices(dofs_per_cell);
std::vector<double> old_u_q(n_q_points);
std::vector<double> old_v_q(n_q_points);
for(const auto &cell: dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
cell_matrix = 0;
cell_rhs = 0;
cell_rhs_1 = 0;
cell_rhs_2 = 0;
fe_values.get_function_values(old_solution_u, old_u_q);
fe_values.get_function_values(old_solution_v, old_v_q);
for(unsigned int q_index = 0; q_index < n_q_points; ++q_index)
{
double old_u = old_u_q[q_index];
double old_v = old_v_q[q_index];
for(unsigned int i = 0; i < dofs_per_cell; ++i)
{
for(unsigned int j = 0; j < dofs_per_cell; ++j)
{
cell_matrix(i, j) += (fe_values.shape_value(i,
q_index) *
fe_values.shape_value(j,
q_index)
+
time_step*time_step*
theta * theta *
fe_values.shape_grad(i,
q_index) *
fe_values.shape_grad(j,
q_index)) *
fe_values.JxW(q_index);
cell_rhs_1(i) +=
/*(M - k^2*theta*(1-theta)*A)U^(n-1)*/
(fe_values.shape_value(i, q_index) *
fe_values.shape_value(j, q_index)
-
time_step*time_step*
theta * (1 - theta)*
fe_values.shape_grad(i, q_index) *
fe_values.shape_grad(j, q_index)
) *
old_u;
cell_rhs_2(i) +=
/*k*M*V^(n-1)*/
time_step *
fe_values.shape_value(i, q_index) *
fe_values.shape_value(j, q_index) *
old_v;
}
cell_rhs(i) += (cell_rhs_1(i) + cell_rhs_2(i) *
fe_values.JxW (q_index));
}
}
cell->get_dof_indices (local_dof_indices);
for(unsigned int i = 0; i < dofs_per_cell; i++)
{
for (unsigned int j = 0; j < dofs_per_cell; j++)
{
matrix_u.add (local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i,j));
}
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
}
{
BoundaryValuesU<dim> boundary_values_u_function;
boundary_values_u_function.set_time (time);
std::map<types::global_dof_index,double> boundary_values;
VectorTools::interpolate_boundary_values (dof_handler,
0,
boundary_values_u_function,
boundary_values);
MatrixTools::apply_boundary_values (boundary_values,
matrix_u,
solution_u,
system_rhs);
}
solve_u ();
std::vector<double> u_q(n_q_points);
for(const auto &cell: dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
cell_matrix = 0;
cell_rhs = 0;
cell_rhs_1 = 0;
cell_rhs_2 = 0;
fe_values.get_function_values(old_solution_u, old_u_q);
fe_values.get_function_values(old_solution_v, old_v_q);
fe_values.get_function_values(solution_u, u_q);
for(unsigned int q_index = 0; q_index < n_q_points; ++q_index)
{
double old_u = old_u_q[q_index];
double old_v = old_v_q[q_index];
double u = u_q[q_index];
for(unsigned int i = 0; i < dofs_per_cell; ++i)
{
for(unsigned int j = 0; j < dofs_per_cell; ++j)
{
cell_matrix(i, j) += (fe_values.shape_value(i,
q_index) *
fe_values.shape_value(j,
q_index)
) *
fe_values.JxW(q_index);
cell_rhs_1(i) +=
/*(M - k^2*theta*(1-theta)*A)U^(n-1)*/
(
fe_values.shape_value(i, q_index) *
fe_values.shape_value(j, q_index)
*
old_v
-
time_step*
(
theta *
fe_values.shape_grad(i, q_index) *
fe_values.shape_grad(j, q_index) *
u
+
(1-theta)*
fe_values.shape_grad(i, q_index) *
fe_values.shape_grad(j, q_index) *
old_u
)
);
}
cell_rhs(i) += (cell_rhs_1(i) *
fe_values.JxW (q_index));
}
}
cell->get_dof_indices (local_dof_indices);
for(unsigned int i = 0; i < dofs_per_cell; i++)
{
for (unsigned int j = 0; j < dofs_per_cell; j++)
{
matrix_u.add (local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i,j));
}
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
}
//***********************************************************************************************************************
{
BoundaryValuesV<dim> boundary_values_v_function;
boundary_values_v_function.set_time (time);
std::map<types::global_dof_index,double> boundary_values;
VectorTools::interpolate_boundary_values (dof_handler,
0,
boundary_values_v_function,
boundary_values);
MatrixTools::apply_boundary_values (boundary_values,
matrix_v,
solution_v,
system_rhs);
}
solve_v ();
output_results ();
std::cout << " Total energy: "
<< (mass_matrix.matrix_norm_square (solution_v) +
laplace_matrix.matrix_norm_square (solution_u)) / 2
<< std::endl;
old_solution_u = solution_u;
old_solution_v = solution_v;
}
}
}
Dulcimer
On Tuesday, January 23, 2018 at 8:27:37 PM UTC+1, Wolfgang Bangerth wrote:
>
> On 01/23/2018 10:35 AM, Dulcimer0909 wrote:
> >
> > If I do go ahead and replace code, so that it does a cell by cell
> assembly, I
> > am a bit lost on how I would store the old_solution (U^(n-1)) for each
> cell
> > and retrieve it during the assembly for the Rhs.
>
> Dulcimer -- can you elaborate? It's not clear to me how the two are
> related.
> In the thread you comment on, the question is about how to assemble the
> matrices, but these matrices do not actually depend on the previous
> solution.
> In any case, the program of course stores the old solution and so any code
> you
> write will have access to it.
>
> So it's not entirely clear to me what you mean in your question :-(
>
> Best
> W.
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: bang...@colostate.edu
> <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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
