Thank you for being kindness!
I followed your suggestions.... but I couldn't resolve the problem yet.
Or did I make other mistake in another places?
>From a slight modification from STEP-4 tutorial , I have been trying solve
Poisson Equation.
Actually, I have manufactured solution to see my error behavior.
But I am testing whether I can really impose dirihlet boundary condition on
the domain boundary.
I would appreciate if you can read my code.....cuz I couldn't resolve the
problem, and having result like this...
although it is the boundary that I impose value of '5'...but it is set at
center....
<https://lh3.googleusercontent.com/-lwNrlyXe0Xg/WCji2e1WY6I/AAAAAAAAA9c/bIEBg6GaTHARYvX-jvghPknmOMaFmnYtQCLcB/s1600/Screen%2BShot%2B2016-11-13%2Bat%2B12.59.38%2BPM.png>
2016년 11월 13일 일요일 오후 2시 43분 46초 UTC-6, Jean-Paul Pelteret 님의 말:
>
> Yes, well this is going to give strange results because you're setting the
> boundary ID of internal faces (x=0.5 is the centreline of your geometry).
> Before you set any boundary ID's, you should always first check to see that
> you're actually on a boundary. The first few tutorials do go through this,
> but I've provided you with a worked example to demonstrate the best
> practise.
>
> Regards,
> J-P
>
>
--
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.
/* ---------------------------------------------------------------------
*
* Copyright (C) 1999 - 2016 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE at
* the top level of the deal.II distribution.
*
* ---------------------------------------------------------------------
*
* Author: Wolfgang Bangerth, University of Heidelberg, 1999
*/
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_handler.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 <fstream>
#include <iostream>
#include <math.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/base/convergence_table.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/numerics/error_estimator.h>
#define Power(x,y) ( std::pow(x,y))
#define Cos(x) ( std::cos(x) )
#define Sin(x) ( std::sin(x) )
#define Sqrt(x) ( std::sqrt(x) )
#define Pi 3.1415926535897932384626
#define E 2.71828182845904523
using namespace dealii;
template <int dim>
class Step4
{
public:
Step4 (const unsigned int degree);
~Step4 ();
void run ();
private:
void make_grid ();
void setup_system();
void assemble_system ();
void solve ();
void refine_grid ();
void error_evaluation();
void output_results () const;
const unsigned int degree;
Triangulation<dim> triangulation;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
ConstraintMatrix constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
double error = 999999;
Vector<double> solution;
Vector<double> system_rhs;
};
/* ---------------------------------------------------------------------
*
* To test code with Manufactured solution
* u=exp(-x^2-y^2);
*
* ---------------------------------------------------------------------*/
template <int dim>
class Solution : public Function<dim>
{
public:
Solution () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double //generating exact solution : u=exp(-x^2-y^2);
Solution<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
Assert (component < this->n_components,
ExcIndexRange (component, 0, this->n_components));
double x=p[0];
double y=p[1];
double return_value;
double S=1; //Determine Strength of Singularity
double z=std::pow(x,2)+std::pow(y,2);
return_value=exp( (-1)*z/S );
return return_value;
}
/* ---------------------------------------------------------------------
*
* To test code with Manufactured solution
* u=exp(-x^2-y^2);
*
* I am trying to impose value of '5' to see whether I can impose Dirihlet BC properly
* ---------------------------------------------------------------------*/
template <int dim>
class BoundaryValues : public Function<dim>
{
public:
BoundaryValues () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double BoundaryValues<dim>::value (const Point<dim> &p,
const unsigned int /*component*/) const
{
// double x=p[0];
// double y=p[1];
// double S=1;
// double z=std::pow(x,2)+std::pow(y,2);
// double return_value;
// return_value=exp( (-1)*z/S );
double return_value=5;
return return_value;
}
/* ---------------------------------------------------------------------
*
* Right Hand side value from Manufactured solution
*
* ---------------------------------------------------------------------*/
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) const
{
double return_value = 0;
double x=p[0]; double y=p[1];
return_value=-4*Power(E,-Power(x,2) - Power(y,2)) + 4*Power(E,-Power(x,2) - Power(y,2))*Power(x,2) +
4*Power(E,-Power(x,2) - Power(y,2))*Power(y,2);
return_value=return_value*(-1);
return return_value;
}
/* ---------------------------------------------------------------------
* OTher Functions...
*
* ---------------------------------------------------------------------*/
template <int dim>
Step4<dim>::Step4 (const unsigned int degree)
:
degree (degree),
fe (degree),
dof_handler (triangulation)
{}
template <int dim>
Step4<dim>::~Step4 ()
{
dof_handler.clear ();
}
template <int dim>
void Step4<dim>::setup_system ()
{
dof_handler.distribute_dofs (fe);
std::cout << " Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< std::endl;
DynamicSparsityPattern dsp(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern (dof_handler, dsp);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit (sparsity_pattern);
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
constraints.clear ();
DoFTools::make_hanging_node_constraints (dof_handler,
constraints);
VectorTools::interpolate_boundary_values (dof_handler,1,
BoundaryValues<dim>(),
constraints);
constraints.close ();
}
template <int dim>
void Step4<dim>::refine_grid ()
{
triangulation.refine_global(1);
}
template <int dim>
void Step4<dim>::assemble_system ()
{
QGauss<dim> quadrature_formula(degree+2);
const MappingQ<dim> mapping (degree);
FEValues<dim> fe_values (mapping, 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);
Vector<double> cell_rhs (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
const RightHandSide<dim> right_hand_side;
std::vector<double> rhs_values (n_q_points);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
fe_values.reinit (cell);
cell_matrix = 0;
cell_rhs = 0;
right_hand_side.value_list (fe_values.get_quadrature_points(), rhs_values);
for (unsigned int q_index=0; q_index<n_q_points; ++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_grad (i, q_index) *
fe_values.shape_grad (j, q_index) *
fe_values.JxW (q_index));
}
cell_rhs(i) += (fe_values.shape_value (i, q_index) *rhs_values [q_index]* fe_values.JxW (q_index));
}
cell->get_dof_indices (local_dof_indices);
constraints.distribute_local_to_global (cell_matrix,cell_rhs,local_dof_indices,system_matrix,system_rhs);
}
}
template <int dim>
void Step4<dim>::solve ()
{
SolverControl solver_control (3000, 1e-12);
SolverCG<> solver (solver_control);
solver.solve (system_matrix, solution, system_rhs,
PreconditionIdentity());
std::cout << " " << solver_control.last_step()
<< " CG iterations needed to obtain convergence."
<< std::endl;
}
template <int dim>
void Step4<dim>::output_results () const
{
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, "solution");
data_out.build_patches ();
std::ofstream output (dim == 2 ?
"solution-2d.vtk" :
"solution-3d.vtk");
data_out.write_vtk (output);
}
template <int dim>
void Step4<dim>::error_evaluation()
{
Vector<float> difference_per_cell (triangulation.n_active_cells());
const MappingQ<dim> mapping (degree);
VectorTools::integrate_difference (mapping, dof_handler,solution,Solution<dim>(),
difference_per_cell,QGauss<dim>(degree+2),VectorTools::L2_norm);
const double L2_error = difference_per_cell.l2_norm();
std::cout << " L2_error : " << L2_error << std::endl;
error=L2_error;
}
/* ---------------------------------------------------------------------
*
* Grid Generation and Boundary designation
*
* ---------------------------------------------------------------------*/
template <int dim>
void Step4<dim>::run ()
{
unsigned int cycle_control=1;
for (unsigned int cycle=0; cycle<cycle_control; ++cycle)
{
if (cycle == 0)
{
GridGenerator::hyper_cube (triangulation, 0.0, 1.0);
triangulation.refine_global (1);
std::cout << "Number of active cells: "<< triangulation.n_active_cells()<< std::endl;
for (typename Triangulation<dim>::active_cell_iterator
cell=triangulation.begin_active();
cell!=triangulation.end(); ++cell)
{
for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
{
// You need this line to ensure that this doesn't affect anything other than bondary faces!
if (cell->face(f)->at_boundary() == true)
{
// This is an internal faces! You never want to set this to anything other than its default value.
// The about check should make sure that we never could execute this part of the code.
if ( std::abs(cell->face(f)->center()[0]-0.5)< 1e-12 )
{
cell->face(f)->set_boundary_id(10);
}
// Boundary faces
static const double tol = 1e-12;
if ( std::abs(cell->face(f)->center()[0]-0.0)< tol )
{
cell->face(f)->set_boundary_id(1);
}
if ( std::abs(cell->face(f)->center()[0]-1.0)< tol )
{
cell->face(f)->set_boundary_id(1);
}
if ( std::abs(cell->face(f)->center()[1]-0.0)< tol )
{
cell->face(f)->set_boundary_id(1);
}
if ( std::abs(cell->face(f)->center()[1]-1.0)< tol )
{
cell->face(f)->set_boundary_id(1);
}
}
}
}
std::ofstream out ("grid-1.ucd");
GridOut grid_out;
grid_out.write_ucd (triangulation, out);
}
else
refine_grid ();
std::cout << " Number of active cells: "<< triangulation.n_active_cells() << std::endl;
setup_system ();
assemble_system ();
solve ();
error_evaluation();
output_results ();
}
}
int main ()
{
Step4<2> laplace_problem_2d1(1);
laplace_problem_2d1.run ();
return 0;
}
