Hi All,
I'm trying to calculate J_tor= x*solution+solution_max/y
solution is already calculated vector from 3 function (setup_system,
assemble_system, solve)
To do this calculation I added the following three function(setup_Jtor,
assmeble_Jtor, solve_Jtor)
template <int dim>
void Step6<dim>::setup_Jtor ()
{
system_rhs.reinit(dof_handler.n_dofs()); //create global space for b
J_tor.reinit (dof_handler.n_dofs()); //create global space for x
system_matrix_Jtor.reinit (sparsity_pattern); //create global space for
A (in Ax=b)
}
template <int dim>
void Step6<dim>::assemble_Jtor ()
{
const QGauss<dim> quadrature_formula(2);
const RightHandSide<dim> right_hand_side;
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();
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell); //create
local space for A
Vector<double> cell_rhs (dofs_per_cell); //create
local space for x
std::vector<double> sol_tmp(n_q_points); //create local space for
solution
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
//by Han(To get psi0)
Vector<double>::iterator max;
max=std::max_element(solution.begin(),solution.end());
for (; cell!=endc; ++cell)
{
fe_values.reinit (cell);
fe_values.get_function_values(solution,sol_tmp);
cell_matrix = 0;
cell_rhs = 0;
for (unsigned int q_index=0; q_index<n_q_points; ++q_index)
{
for(unsigned int i=0; i<dofs_per_cell; ++i) //local calculation for A
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));
for(unsigned int i=0; i<dofs_per_cell; ++i) //local calculation for b
cell_rhs(i) += (fe_values.shape_value(i, q_index)*
right_hand_side.value(fe_values.quadrature_point(q_index),
sol_tmp[q_index],*max) *
fe_values.JxW (q_index));
}
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i) //move calculated
values from local to global
for (unsigned int j=0; j<dofs_per_cell; ++j)
system_matrix_Jtor.add (local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i,j));
for (unsigned int i=0; i<dofs_per_cell; ++i)
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
std::map<types::global_dof_index,double> boundary_values;
VectorTools::interpolate_boundary_values (dof_handler,
0,
ZeroFunction<dim>(),
boundary_values);
MatrixTools::apply_boundary_values (boundary_values,
system_matrix_Jtor,
solution,
system_rhs);
}
template <int dim>
void Step6<dim>::solve_Jtor ()
{
SolverControl solver_control (2000, 1e-12);
SolverCG<> solver (solver_control);
PreconditionSSOR<> preconditioner;
preconditioner.initialize(system_matrix_Jtor, 1.2);
solver.solve (system_matrix_Jtor, J_tor, system_rhs,
preconditioner);
}
In other words, I just calculated (pi_i,pi_j) x J_tor_j=(rhs_i,pi_i)
(where pi is test function)
It worked, the results are as I expected. But someone said to me it is very
inefficient way.
I agree with that opinion
But I have to access to Point(x and y), solution and maximum value of
solution to calculate J_tor
And, I also have to visualize J_tor using Visit.
Could you let me know how can I solve J_tor more efficiently using the
above three information?
Thank you.
Kyusik.
--
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.
