And this is my code:
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/function.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/convergence_table.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/tensor_function.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/solution_transfer.h>
#include <math.h>
#include <iostream>
#include <fstream>
#include <sstream>
#define pi 3.1415926
namespace Case1
{
using namespace dealii;
template <int dim>
class PhaseEquation
{
public:
PhaseEquation ();
void run ();
private:
void setup_system ();
void solve_phi_star ();
void assemble_system ();
// void output_results () const;
//double compute_errors () const; //i need to look for this example
Triangulation<dim> triangulation;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
ConstraintMatrix constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> mass_matrix;
SparseMatrix<double> laplace_matrix;
SparseMatrix<double> matrix_phi;
SparseMatrix<double> mass_bar_matrix;
Vector<double> solution_phi, solution_phi_star;
Vector<double> old_solution_phi;
Vector<double> system_rhs;
double time_step, time;
const double eps;//, s;
unsigned int timestep_number;
// const double eps;
};
template <int dim>
PhaseEquation<dim>::PhaseEquation () :
fe (1),
dof_handler (triangulation),
time_step (1./64),
time (time_step),
timestep_number (1),
eps (0.03)
{}
template <int dim>
class Solution : public Function<dim>
// protected SolutionBase<dim>
{
public:
Solution () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
virtual Tensor<1,dim> gradient (const Point<dim> &p,
const unsigned int component = 0)
const;
};
// The actual definition of the values and gradients of the exact solution
template <int dim>
double Solution<dim>::value (const Point<dim> &p,
const unsigned int) const
{
double return_value = 0;
/*for (unsigned int i=0; i<this->n_source_centers; ++i)
{
const Tensor<1,dim> x_minus_xi = p - this->source_centers[i];
return_value += std::exp(-x_minus_xi.norm_square() /
(this->width * this->width));
}*/
/*return 0.5 * (1 - std::tanh((p[0] - (3 / (std::sqrt(2) * eps) *
this->get_time())) / (2 * std::sqrt(2) * eps)))*/
return return_value;
}
template <int dim>
Tensor<1,dim> Solution<dim>::gradient (const Point<dim> &p,
const unsigned int) const
{
Tensor<1,dim> return_value;
/*for (unsigned int i=0; i<this->n_source_centers; ++i)
{
const Tensor<1,dim> x_minus_xi = p - this->source_centers[i];
// For the gradient, note that its direction is along (x-x_i), so we
// add up multiples of this distance vector, where the factor is
given
// by the exponentials.
return_value += (-2 / (this->width * this->width) *
std::exp(-x_minus_xi.norm_square() /
(this->width * this->width)) *
x_minus_xi);
}*/
return_value[0] = 0;
return_value[1] = 0;
return return_value;
}
template <int dim>
class InitialValue : public Function<dim>
{
public:
InitialValue () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double InitialValue<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
(void) component;
Assert(component == 0, ExcIndexRange(component, 0, 1));
const double eps = 0.03;
return 0.5 * (1 - std::tanh((p[0] - (3 / (std::sqrt(2) * eps) * 0)) /
(2 * std::sqrt(2) * eps)));
}
// Secondly, we have the right hand side forcing term. Boring as we are,
we
// choose zero here as well:
template <int dim>
class RightHandSide : public Function<dim>
{
public:
RightHandSide () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double RightHandSide<dim>::value (const Point<dim> &/*p*/,
const unsigned int component) const
{
(void) component;
Assert(component == 0, ExcIndexRange(component, 0, 1));
return 0;
}
// Finally, we have boundary values for $u$ and $v$. They are as described
// in the introduction, one being the time derivative of the other:
template <int dim>
class BoundaryValue : public Function<dim>
{
public:
BoundaryValue () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double BoundaryValue<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
(void) component;
Assert(component == 0, ExcIndexRange(component, 0, 1));
const double eps = 0.03;
if ((this->get_time() <= 1) &&
((p[0] - (-1)) <= 1e-12) &&
((p[1] - (-1)) <= 1e-12) )
return 0.5 * (1 - std::tanh((-1 - (3 / (std::sqrt(2) * eps) *
this->get_time())) / (2 * std::sqrt(2) * eps)));//std::sin
(this->get_time() * 4 * numbers::PI);
//else
//return 0;
}
template <int dim>
void PhaseEquation<dim>::setup_system ()
{
GridGenerator::hyper_cube (triangulation, -1, 1);
triangulation.refine_global (7);
typename Triangulation<dim>::cell_iterator
cell = triangulation.begin (),
endc = triangulation.end();
for (; cell!=endc; ++cell){
for (unsigned int face_number=0;
face_number<GeometryInfo<dim>::faces_per_cell;
++face_number){
if (((std::fabs(cell->face(face_number)->center()(0) - (-1)) >
1e-12)
||
(std::fabs(cell->face(face_number)->center()(1) - (-1)) >
1e-12)))
/*||
(()
&&
()
))*/
cell->face(face_number)->set_boundary_id (1);}
}
std::cout << "Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl;
dof_handler.distribute_dofs (fe);
std::cout << "Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< std::endl
<< std::endl;
DynamicSparsityPattern dsp(dof_handler.n_dofs());//,
dof_handler.n_dofs());
DoFTools::make_sparsity_pattern (dof_handler, dsp);
sparsity_pattern.copy_from (dsp);
// Then comes a block where we have to initialize the 3 matrices we need
mass_matrix.reinit (sparsity_pattern);
laplace_matrix.reinit (sparsity_pattern);
matrix_phi.reinit (sparsity_pattern);
mass_bar_matrix.reinit (sparsity_pattern);
//MatrixCreator::create_mass_matrix (dof_handler, QGauss<dim>(3),
// mass_matrix);
//MatrixCreator::create_laplace_matrix (dof_handler, QGauss<dim>(3),
// laplace_matrix);
// The rest of the function is spent on setting vector sizes to the
solution_phi.reinit (dof_handler.n_dofs());
solution_phi_star.reinit (dof_handler.n_dofs());
old_solution_phi.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
constraints.close ();
}
template <int dim>
void PhaseEquation<dim>::assemble_system ()
{
QGauss<dim> quadrature_formula(3);
QGauss<dim-1> face_quadrature_formula(3);
const unsigned int n_q_points = quadrature_formula.size();
const unsigned int n_face_q_points = face_quadrature_formula.size();
const unsigned int dofs_per_cell = fe.dofs_per_cell;
FullMatrix<double> cell_mass_matrix (dofs_per_cell, dofs_per_cell);
FullMatrix<double> cell_laplace_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula,
update_values |
update_quadrature_points |
update_normal_vectors |
update_JxW_values);
const RightHandSide<dim> right_hand_side;
std::vector<double> rhs_values (n_q_points);
const Solution<dim> exact_solution;
// Now for the main loop over all cells. This is mostly unchanged from
// previous examples, so we only comment on the things that have
changed.
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
cell_mass_matrix = 0;
cell_laplace_matrix = 0;
cell_rhs = 0;
fe_values.reinit (cell);
right_hand_side.value_list (fe_values.get_quadrature_points(),
rhs_values);
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
// The first thing that has changed is the bilinear form. It
// now contains the additional term from the Helmholtz
// equation:
cell_mass_matrix(i,j) += (fe_values.shape_value(i,q_point) *
fe_values.shape_value(j,q_point)) *
fe_values.JxW(q_point);
cell_laplace_matrix(i,j) +=(fe_values.shape_grad(i,q_point) *
fe_values.shape_grad(j,q_point)) *
fe_values.JxW(q_point);
}
cell_rhs(i) += (fe_values.shape_value(i,q_point) *
rhs_values [q_point] *
fe_values.JxW(q_point));
}
for (unsigned int face_number=0;
face_number<GeometryInfo<dim>::faces_per_cell; ++face_number)
{
if (cell->face(face_number)->at_boundary()
&&
(cell->face(face_number)->boundary_id() == 1))
{
// If we came into here, then we have found an external face
// belonging to Gamma2. Next, we have to compute the values of
// the shape functions and the other quantities which we will
// need for the computation of the contour integral. This is
// done using the <code>reinit</code> function which we
already
// know from the FEValue class:
fe_face_values.reinit (cell, face_number);
// And we can then perform the integration by using a loop
over
// all quadrature points.
//
// On each quadrature point, we first compute the value of the
// normal derivative. We do so using the gradient of the exact
// solution and the normal vector to the face at the present
// quadrature point obtained from the
// <code>fe_face_values</code> object. This is then used to
// compute the additional contribution of this face to the
right
// hand side:
for (unsigned int q_point=0; q_point<n_face_q_points;
++q_point)
{
const double neumann_value
= (exact_solution.gradient
(fe_face_values.quadrature_point(q_point)) *
fe_face_values.normal_vector(q_point));
for (unsigned int i=0; i<dofs_per_cell; ++i){
cell_rhs(i) += (neumann_value *
fe_face_values.shape_value(i,q_point) *
fe_face_values.JxW(q_point));}
}
}
}
// Now that we have the contributions of the present cell, we can
// transfer it to the global matrix and right hand side vector, as
in
// the examples before:
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){
mass_matrix.add (local_dof_indices[i],
local_dof_indices[j],
cell_mass_matrix(i,j));
laplace_matrix.add (local_dof_indices[i],
local_dof_indices[j],
cell_laplace_matrix(i,j));
system_rhs(local_dof_indices[i]) += cell_rhs(i);}
}
mass_bar_matrix = 0;
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
mass_bar_matrix.add (local_dof_indices[i],
local_dof_indices[i],
mass_matrix(i,j));
}
}
//hanging_node_constraints.condense (system_matrix);
//hanging_node_constraints.condense (system_rhs);
}
}
// without:
template <int dim>
void PhaseEquation<dim>::solve_phi_star ()
{
SolverControl solver_control (1000,
1e-8*system_rhs.l2_norm());
SolverCG<> cg (solver_control);
cg.solve (mass_bar_matrix, solution_phi_star, system_rhs,
PreconditionIdentity());
std::cout << " u-equation: " << solver_control.last_step()
<< " CG iterations."
<< std::endl;
}
/* template <int dim>
void WaveEquation<dim>::compute_errors () const
{
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution_u, "U");
data_out.add_data_vector (solution_v, "V");
data_out.build_patches ();
const std::string filename = "solution-" +
Utilities::int_to_string (timestep_number,
3) +
".gnuplot";
std::ofstream output (filename.c_str());
data_out.write_gnuplot (output);
}
*/
template <int dim>
void PhaseEquation<dim>::run ()
{
setup_system();
assemble_system ();
VectorTools::project (dof_handler, constraints, QGauss<dim>(3),
InitialValue<dim>(),
old_solution_phi);
Vector<double> tmp (solution_phi.size());
//Vector<double> forcing_terms (solution_phi.size());
for (; time<=5; time+=time_step, ++timestep_number)
{
std::cout << "Time step " << timestep_number
<< " at t=" << time
<< std::endl;
mass_matrix.vmult (system_rhs, old_solution_phi);
{
BoundaryValue<dim> boundary_value_phi_function;
boundary_value_phi_function.set_time (time);
std::map<types::global_dof_index,double> boundary_value;
VectorTools::interpolate_boundary_values (dof_handler,
0,
boundary_value_phi_function,
boundary_value);
matrix_phi.copy_from (mass_bar_matrix);
matrix_phi.add (time_step, laplace_matrix);
MatrixTools::apply_boundary_values (boundary_value,
matrix_phi,
solution_phi_star, // star
is this? bc is not for star
system_rhs);
solve_phi_star ();
solution_phi[0] = solution_phi_star[0] * (1 / (std::exp(((-2) *
time_step) / (eps * eps)) + (solution_phi_star[0] * solution_phi_star[0] +
solution_phi_star[1] * solution_phi_star[1]) * (1 - std::exp(((-2) *
time_step) / (eps * eps)))));//
solution_phi[1] = solution_phi_star[1] * (1 / (std::exp(((-2) *
time_step) / (eps * eps)) + (solution_phi_star[0] * solution_phi_star[0] +
solution_phi_star[1] * solution_phi_star[1]) * (1 - std::exp(((-2) *
time_step) / (eps * eps)))));
//output_results ();
//compute_errors ();
old_solution_phi = solution_phi;
}
}
////compute_errors ();
}
int main ()
{
//try
//{
using namespace dealii;
using namespace Case1;
PhaseEquation<2> phase_equation_solver;
phase_equation_solver.run ();
// }
/*catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
*/
return 0;
}
}
I think that my code has something wrong, because the examples of deal.ii
can be run, maybe it is not the machine's problem. But I don't know what's
wrong in my code?
chucui
在 2018年5月27日星期日 UTC+8下午2:41:00，chucu...@gmail.com写道：
>
> Hi, everyone,
>
> I am trying to implement some easier case by using deal.ii, and I have a
> code and debug it, after cmake and when make it, the results on the
> Terminal are like this:
>
>
> <https://lh3.googleusercontent.com/-iVRlUfiS2A0/WwpRJlFrGxI/AAAAAAAAABE/EOWp447IRAYWBGpr3miMALUwEJyu1jzgwCLcBGAs/s1600/1.JPG>
> I type: yum list | grep libgfortran
> the result is:
>
>
> <https://lh3.googleusercontent.com/-fbmkFDbe3ds/WwpR0a6in1I/AAAAAAAAABM/kpmB_FK4iCQTgEhfBmc9Qh1j8Oi59APSQCLcBGAs/s1600/2.JPG>
> I don't know what happen in my machine, but when I run the examples of
> deal.ii, every thing is OK.
>
> What should I do next, and how can I debug it successfully? Thank you too
> much!
>
> chucui
>
--
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.
