Dear Jonathan,
Thank you very much! Firstly, My vtk file's names are allowed yo have more
than three digits. My code of writting vtk file is:
template <int dim>
void MyProblem<dim>::output_results_U (const unsigned int
timestep_number) const
{
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler_U);
data_out.add_data_vector (old_U, "U");
data_out.build_patches (degree+3);
std::ostringstream filename;
filename << "U-"
<< Utilities::int_to_string(timestep_number)
<< ".vtk";
std::ofstream output (filename.str().c_str());
data_out.write_vtk (output);
}
And I have got VTK files such as 'U-6400.vtk'
I'm sorry for my discription without codes. For example, the system I want
to solve is just like:
[image: ol.JPG]
Then the weak form in each time step is:
[image: 哦了.JPG]
The code to assemble system matrix and right hand side is:
template <int dim>
void MyProblem<dim>::assemble_system_U (const Vector<double> old_U)
{
system_matrix_U=0;
system_rhs_U = 0;
QGauss<dim> quadrature_formula(degree+2);
FEValues<dim> fe_U_values (fe_U, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_hessians |
update_JxW_values);
const unsigned int dofs_per_cell_U = fe_U.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
Vector<double> local_rhs (dofs_per_cell_U);
FullMatrix<double> local_matrix (dofs_per_cell_U, dofs_per_cell_U);
std::vector<types::global_dof_index> local_dof_indices
(dofs_per_cell_U);
std::vector<double> old_U_values(n_q_points);
typename DoFHandler<dim>::active_cell_iterator
cell_U = dof_handler_U.begin_active(),
endc_U = dof_handler_U.end();
for (; cell_U!=endc_U; ++cell_U)
{
local_rhs = 0;
local_matrix=0;
fe_U_values.reinit (cell_U);
fe_U_values.get_function_values (old_U, old_U_values);
for (unsigned int q=0; q<n_q_points; ++q)
{
const double old_U = old_U_values[q];
for (unsigned int i=0; i<dofs_per_cell_U; ++i)
{
const double psi_i_U = fe_U_values.shape_value (i, q);
for (unsigned int j=0; j<dofs_per_cell_U; ++j)
{
const double psi_j_U = fe_U_values.shape_value (j,
q);
local_matrix(i,j) += (psi_j_U * psi_i_U * (1./time_step))
* fe_U_values.JxW(q);
}
local_rhs(i) += (old_U * psi_i_U * (1./time_step)
+ F_Function(old_U) * psi_i_U ) *
fe_U_values.JxW(q);
}
}
cell_U->get_dof_indices (local_dof_indices);
constraints_U.distribute_local_to_global(local_matrix, local_rhs,
local_dof_indices, system_matrix_U, system_rhs_U);
}
}
Then solve new_U i.e. U_{n+1}, and update the U vector.
template <int dim>
void MyProblem<dim>::solve_U ()
{
SparseDirectUMFPACK A_direct;
A_direct.initialize(system_matrix_U);
A_direct.vmult (new_U, system_rhs_U);
}
template <int dim>
void MyProblem<dim>::run ()
{
[...]
VectorTools::project (dof_handler_U, constraints_U,
QGauss<dim>(degree+2),
InitialValues_U<dim>(),
old_U);
int timestep_number = 1;
for (; time<=100000000.0; time+=time_step, ++timestep_number)
{
assemble_system_U (old_U);
solve_U ();
constraints_U.distribute (new_U);
old_U = new_U;
output_results_U (timestep_number);}}
But if the code is broken suddenly ,at time step n=1000 for example, how to
restart it from time step t_n = 1001, not t_n=1? And in each time step, I
have a output VTK file of solution U (U-n.vtk), but in the code I need
‘const Vector<double> old_U’ in fact. How to convert a U-1000.vtk file into
a 'const Vector<double> old_U', then use it for 'assemble_system_U (const
Vector<double> old_U)'? Thank you very much!
Best,
Chucui
在 2018年11月5日星期一 UTC+8下午1:25:58，[email protected]写道：
>
> If you're using the timestep number in the file names of your vtk files,
> are you formatting the filename to allow more than three digits? (e.g.
> solution-999.vtk versus solution-0999.vtk.)
>
> It's hard to tell without seeing any of your code.
> Cheers,
> Jonathan
>
>
--
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.
