2020-12-10 수정
2021-03-10 수정
tutorial3의 기본적인 구조는 아래와 같다.
이 프로그램은 유한 요소 방법을 사용하여 해결할 수있는 가장 간단한 방정식이지만,
대부분의 유한 요소 프로그램의 기본 구조는 다음과 같다.
거의 모든 본질적으로 따를 템플릿 역할을합니다.
class Step3
{
public:
Step3 ();
void run ();
private:
void make_grid ();
void setup_system ();
void assemble_system ();
void solve ();
void output_results () const;
Triangulation<2> triangulation;
FE_Q<2> fe;
DoFHandler<2> dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};이것은 데이터 캡슐화 data encapsulation개체 지향 프로그래밍은
외부에 액세스 할 수없는 개인 멤버에 클래스의 거의 모든
내부 세부 사항을 숨기려고한다.
개체 지향 프로그래밍(OOP)에서 캡슐화는 해당 데이터에서 작동하는 메서드와 데이터 번들링 또는 개
체의 일부 구성 요소에 대한 직접 액세스 제한을 말합니다. ==> 말 그대로 숨겨서 표현하고 싶다.
멤버 변수
객체 지향 프로그래밍에서 멤버 변수(member variable) 또는 멤버 필드는 특정 객체와 연결된 변수의 하나이며, 해당 변수의 모든 메소드(멤버 함수)에 접근이 가능하다. 클래스 기반 언어에서 이들은 두 종류로 구별된다: 모든 인스턴스의 클래스와 공유되는 변수의 사본이 하나만 있을 경우 이를 클래스 변수나 정적 멤버 변수로 부른다. 클래스의 각 인스턴스가 자신만의 변수 복사본을 소유하고 있는 경우 해당 변수는 인스턴스 변수라 부른다.[1]
객체 지향 프로그래밍에서 멤버 변수(member variable) 또는 멤버 필드는 특정 객체와 연결된 변수의 하나이며, 해당 변수의 모든 메소드(멤버 함수)에 접근이 가능하다. 클래스 기반 언어에서 이들은 두 종류로 구별된다: 모든 인스턴스의 클래스와 공유되는 변수의 사본이 하나만 있을 경우 이를 클래스 변수나 정적 멤버 변수로 부른다. 클래스의 각 인스턴스가 자신만의 변수 복사본을 소유하고 있는 경우 해당 변수는 인스턴스 변수라 부른다.[1]
객체 지향 프로그래밍에서 멤버 변수(member variable) 또는 멤버 필드는 특정 객체와 연결된 변수의 하나이며, 해당 변수의 모든 메소드(멤버 함수)에 접근이 가능하다. 클래스 기반 언어에서 이들은 두 종류로 구별된다: 모든 인스턴스의 클래스와 공유되는 변수의 사본이 하나만 있을 경우 이를 클래스 변수나 정적 멤버 변수로 부른다. 클래스의 각 인스턴스가 자신만의 변수 복사본을 소유하고 있는 경우 해당 변수는 인스턴스 변수라 부른다.[1]
객체 지향 프로그래밍에서 멤버 변수(member variable) 또는 멤버 필드는 특정 객체와 연결된 변수의 하나이며, 해당 변수의 모든 메소드(멤버 함수)에 접근이 가능하다. 클래스 기반 언어에서 이들은 두 종류로 구별된다: 모든 인스턴스의 클래스와 공유되는 변수의 사본이 하나만 있을 경우 이를 클래스 변수나 정적 멤버 변수로 부른다. 클래스의 각 인스턴스가 자신만의 변수 복사본을 소유하고 있는 경우 해당 변수는 인스턴스 변수라 부른다.[1]
#include <iostream>
class Foo {
int bar; // 멤버 변수
public:
void setBar (const int newBar) { bar = newBar; }
};
int main () {
Foo rect; // 지역 변수
return 0;
}
멤버 변수로 시작해 봅시다:
Triangulation<2> triangulation;
FE_Q<2> fe;
DoFHandler<2> dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;글머리 기호에 위에서 설명한 빌딩 블록,
즉 Triangulation 과 DoFHandler object, 그리고 사용하려는 셰이프 함수의 종류를
설명하는 유한 요소 오브젝트를 따릅니다.
두 번째 개체 그룹은 선형 대수와 관련이 있습니다: 시스템 매트릭스및 오른쪽뿐만 아니라 솔루션 벡터, 그리고 행렬의 희소성 패턴을 설명하는 객체. 이것은 이 모든 클래스의 요구 사항(고정된 PDE에 대한 모든 솔버가 요구하는 필수 요소)이며 전체 프로그램 전체에서 생존해야 합니다.
이와 는 대조적으로 Assembly에 필요한
FEValue object는 어셈블리 전체에만 필요하므로 이를 수행하는 함수의 로컬 개체로 만들고 마지막에 다시 파괴합니다.
둘째, 멤버 함수를 살펴보겠습니다.
public:
Step3 ();
void run ();
private:
void make_grid ();
void setup_system ();
void assemble_system ();
void solve ();
void output_results () const;이들은 모든 튜토리얼 프로그램이 사용하는 공통 구조를 형성합니다.
- make_grid(): '전처리 함수'라고 부를 수 있는 것입니다. 이름에서 알 수 있듯이 Triangulation을 저장하는 개체를 설정합니다. 이후 예제에서는 경계 조건, 형상 등을 처리할 수도 있습니다.
- setup_system(): 다음 문제를 해결하는 데 필요한 다른 모든 데이터 구조가 설정된 기능입니다.
특히, DoFHandler 오브젝트를 초기화하고 선형 대수와 관련이 있는 다양한 개체의 크기를 올바르게 조정합니다. 이 함수는 시간 종속 프로그램에서 메시가 적절하게 정제될(mesh is adaptively refined) 때마다 적어도 몇 시간마다 호출될 수 있기 때문에 위의 전처리 기능과 분리되는 경우가 많습니다(스텝-6에서수행하는 방법을 볼 수 있음). 한편, 위의 전처리 함수에서 메시 자체를 설정하는 것은 프로그램의 시작 부분에서 한 번만 수행되므로 자체 함수로 분리됩니다. ==> DoFHandler는 메시가 바뀔 때 마다 호출될 수 있어서 위의 전처리 기능과 분리 되는 경우가 많음 - assemble_system(): 행렬과 오른쪽의 내용이 계산되는 곳입니다. 이 선형 시스템으로 작업을 수행하는 것은 항목을 계산하는 것과 개념적으로 매우 다르기 때문에 다음 함수(solve)와 분리합니다.
- solve(): U 솔루션을 계산한다. 선형 시스템의 AU=F. 현재 프로그램에서는 행렬이 너무 간단하기 때문에 간단한 작업이지만 문제가 더 이상 사소한 것이 아닐 때마다 프로그램 크기의 중요한 부분이 될 것입니다 (단계-20, step-22또는 Step-31을 참조)
- output_results(): 마지막으로 솔루션을 계산하면 솔루션을 사용하여 뭔가를 하고 싶을 것입니다. 예를 들어 시각화할 수 있는 형식으로 출력하거나 관심 있는 수량을 계산할 수 있습니다: 열 교환기의 열 플럭스, 날개의 공기 마찰 계수, 최대 교량 하중 또는 단순히 지점에서 수치 솔루션의 값. 따라서 이 함수는 솔루션을 후처리할 수 있는 장소입니다.
이 모든 것은 단일 공용 함수(생성자 제외)에 의해 함께 유지됩니다(즉, 함수).
이 형식의 개체가 생성되는 위치에서 호출되는 것이며, 다른 모든 함수를 적절한 순서로 호출하는 것입니다.
이 작업을 다른 모든 함수를 호출하는 대신 함수에 캡슐화하면 이 클래스 내에서 우려 분리가 구현되는 방식을 변경할 수 있습니다. 예를 들어 함수 중 하나가 너무 커지면 함수를 두 개로 나눌 수 있으며 결과적으로 변경에 대해 염려해야 하는 유일한 장소는 이 클래스 내에 있으며 다른 곳에서는 그렇지 않습니다.run()run()main()
위에서 언급했듯이 함수 이름 철자가 변형된 이 일반적인 구조는 볼 수 있지만, 본질적으로 이러한 기능 분리 순서는 다음 자습서 프로그램의 많은 부분에서 다시 볼 수 있습니다.
이 모든 것은 단일 공용 함수(생성자 제외)에 의해 함께 유지됩니다(즉, 함수).
이 형식의 개체가 생성되는 위치에서 호출되는 것이며, 다른 모든 함수를 적절한 순서로 호출하는 것입니다.
이 작업을 다른 모든 함수를 호출하는 대신 함수에 캡슐화하면 이 클래스 내에서 우려 분리가 구현되는 방식을 변경할 수 있습니다. 예를 들어 함수 중 하나가 너무 커지면 함수를 두 개로 나눌 수 있으며 결과적으로 변경에 대해 염려해야 하는 유일한 장소는 이 클래스 내에 있으며 다른 곳에서는 그렇지 않습니다.run()run()main()
위에서 언급했듯이 함수 이름 철자가 변형된 이 일반적인 구조는 볼 수 있지만, 본질적으로 이러한 기능 분리 순서는 다음 자습서 프로그램의 많은 부분에서 다시 볼 수 있습니다.
이 모든 것은 단일 공용 함수(생성자 제외)에 의해 함께 유지됩니다(즉, 함수).
이 형식의 개체가 생성되는 위치에서 호출되는 것이며, 다른 모든 함수를 적절한 순서로 호출하는 것입니다.
이 작업을 다른 모든 함수를 호출하는 대신 함수에 캡슐화하면 이 클래스 내에서 우려 분리가 구현되는 방식을 변경할 수 있습니다. 예를 들어 함수 중 하나가 너무 커지면 함수를 두 개로 나눌 수 있으며 결과적으로 변경에 대해 염려해야 하는 유일한 장소는 이 클래스 내에 있으며 다른 곳에서는 그렇지 않습니다.run()run()main()
위에서 언급했듯이 함수 이름 철자가 변형된 이 일반적인 구조는 볼 수 있지만, 본질적으로 이러한 기능 분리 순서는 다음 자습서 프로그램의 많은 부분에서 다시 볼 수 있습니다.
이 모든 것은 단일 공용 함수(생성자 제외)에 의해 함께 유지됩니다(즉, 함수).
이 형식의 개체가 생성되는 위치에서 호출되는 것이며, 다른 모든 함수를 적절한 순서로 호출하는 것입니다.
이 작업을 다른 모든 함수를 호출하는 대신 함수에 캡슐화하면 이 클래스 내에서 우려 분리가 구현되는 방식을 변경할 수 있습니다. 예를 들어 함수 중 하나가 너무 커지면 함수를 두 개로 나눌 수 있으며 결과적으로 변경에 대해 염려해야 하는 유일한 장소는 이 클래스 내에 있으며 다른 곳에서는 그렇지 않습니다.run()run()main()
위에서 언급했듯이 함수 이름 철자가 변형된 이 일반적인 구조는 볼 수 있지만, 본질적으로 이러한 기능 분리 순서는 다음 자습서 프로그램의 많은 부분에서 다시 볼 수 있습니다.
라이브러리 임포트
Triangulation과 DoF선언
#include <deal.II/grid/tria.h>
#include <deal.II/dofs/dof_handler.h>
그리드를 만드는 함수가 선언되는 파일
#include <deal.II/grid/grid_generator.h>
모든 셀의 루프에 필요한 클래스와 셀 오브젝트에서 정보를 가져오는 클래스가 포함되어 있습니다.
처음 두 cell에서 기하학적 정보를 얻기 위해 전에 사용 되었습니다.
마지막 하나는 새로운 및 셀에 로컬 자유의 정도에 대 한 정보를 제공 합니다.
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_accessor.h>
Lagrange 보간 유한 요소에 대한 설명이 포함
#include <deal.II/fe/fe_q.h>
희소성 행렬의 희소성 패턴을 만드는 데 필요
#include <deal.II/dofs/dof_tools.h>
아래 코드는 matrix를 Assembly하는데 사용된다.
각셀에서 quadrature를 사용한다. (quadrature is a historical term which means the process of determining area.)
#include <deal.II/fe/fe_values.h>
#include <deal.II/base/quadrature_lib.h>
경계 값 처리를 위해 필요한 파일
#include <deal.II/base/function.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
두 번째로 마지막 그룹인 포함 파일은 Laplace 방정식의 유한 요소 이산으로부터 발생하는 방정식 시스템을 해결하기 위해 사용하는 선형 대수입니다. 각 셀에서 방정식 시스템을 로컬로 조립하고 결과를 희소 매트릭스로 전송하기 위해 벡터와 전체 행렬을 사용합니다. 그런 다음 Conjugate Gradient solver를 사용하여 사전 조건이 필요한 문제를 해결합니다(이 프로그램에서는 아무 것도 하지 않는 identity preconditioner를 사용하지만 어쨌든 파일을 포함해야 합니다).
#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>
Finally, this is for output to a file and to the console:
#include <deal.II/numerics/data_out.h>
#include <fstream>
#include <iostream>
글로벌 범위로 deal.II 네임스페이스를 가져오는 것
using namespace dealii;
The Step3 class
클래스의 public 부분은 다소 짧습니다.
생성자와 외부에서 호출되고 main 함수와 같은 역할을하는 함수 run 있습니다.
이 클래스의 어떤 작업이 어떤 순서로 실행될 것인지 조정합니다.
클래스의 다른 모든 것, 즉 실제로 어떤 일을하는 모든 함수는 클래스의 private 섹션에 있습니다.
class Step3
{
public:
Step3();
void run();
대부분의 name이 suggest한 멤버 함수가 있습니다.
외부에서 호출 할 필요가 없기 때문에 이 클래스에 private로 설정됩니다.
private:
void make_grid();
void setup_system();
void assemble_system();
void solve();
void output_results() const;
삼각 분할과 자유도의 전역 번호 매기기를 설명하는 멤버 변수가 있습니다
(이 클래스의 생성자에서 유한 요소의 정확한 다항식 정도를 지정합니다).
Triangulation<2> triangulation;
FE_Q<2> fe;
DoFHandler<2> dof_handler;
Laplace 방정식의 discretization으로 인한 시스템 매트릭스의 희소성 패턴 및 값에 대한 변수
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
오른쪽 및 솔루션 벡터를 보유하는 변수
Vector<double> solution;
Vector<double> system_rhs;
};
Step3::Step3
여기에 생성자(construct)가 온다. 다각도를 나타내는 유한 요소 개체에 대한 매개 변수에 의해 표시됨)을 지정하고
dof_handler 변수를 사용하는 Triangulation에 연결하는 것은 먼저 그다지 많지 않습니다.
(Triangulation은 현재 메쉬로 설정되지 않지만 DoFHandler는 상관하지 않습니다 : 그것은 단지 어떤 Triangulation과 연관될지 알고 싶어하며, distribute_dofs() 함수를 사용하여 메시에 자유의 정도를 배포하려고 하면 실제 메시에 만 관심을 기울이기 시작합니다.)
Step3 클래스의 다른 모든 멤버 변수에는 원하는 모든 작업을 수행하는 기본 생성자가 있습니다.
Step3::Step3()
: fe(1)
, dof_handler(triangulation)
{}
Step3::make_grid
이제 가장 먼저 해야 할 일은 계산을 하고 각 정점을 자유의 수준으로 번호 매기고자 하는 Triangulation을 생성하는 것입니다. 우리는 각각 전에 단계 1및 단계 2에서 이 두 단계를 보았습니다.
이 함수는 첫 번째 부분을 수행하여 메시를 만듭니다.
우리는 그리드를 만들고 모든 셀을 5 번 다듬습니다.
초기 그리드 이후(정사각형입니다. [−1,1]×[-1,1])는 하나의 셀로 만 구성되며,
최종 그리드는 32배 32셀을 가지며, 총 1024개의 셀이다.
1024가 올바른 숫자인지 확실하지 않은가요?
우리는 Triangulation에 기능을 사용하여 세포의 수를 아웃하여 확인할 수 있습니다.n_active_cells()
void Step3::make_grid()
{
GridGenerator::hyper_cube(triangulation, -1, 1);
triangulation.refine_global(5);
std::cout << "Number of active cells: " << triangulation.n_active_cells()
<< std::endl;
}
NOTE
Triangulation::n_cells()대신 Triangulation::n_active_cells() 함수라고 부른다.
여기서 active은 더 이상 정제되지 않은 cell 의미합니다.
초기 그리드를 구성하는 하나의 셀까지 더 많은 셀,
즉 가장 미세한 셀의 상위 셀, 상위 셀 등이 있기 때문에 형용사 "active"를 강조합니다.
물론, 다음 거친 수준에서, 세포의 수는 최고 수준에 세포의 1 분기,
즉 256, 다음 64, 16, 4, 및 1. 위의 코드에서 대신 호출하면
결과적으로 1365의 값을 얻게 됩니다.
한편, cell의 개수(활성 셀 수와 는 반대로)는 전형적으로
많은 관심을 받지 않으므로 이를 print할 이유가 없다.triangulation.n_cells()
Step3::setup_system
다음으로 모든 자유도를 열거(enumerate)하고 시스템 데이터를 보관할 행렬과 벡터 객체를 설정합니다.
Enumerating(열거)는 step-2 예제에서 본 것처럼 DoFHandler::distribute_dofs()를 사용하여 수행됩니다.
FE_Q클래스를 사용하고 생성자, 즉 쌍 선형 요소(bilinear elements)에서 다항식 차수를 1로 설정 했으므로
각 정점에 1 개의 자유도를 연관시킵니다.
출력을 생성하는 동안 얼마나 많은 자유도(degrees of freedom)가 생성되는지 살펴 보겠습니다.
void Step3::setup_system()
{
dof_handler.distribute_dofs(fe);
std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
각 정점(vertex)마다 하나의 DoF가 있어야 합니다. 3
2배 32그리드를 가지고 있기 때문에 DoF의 수는 33배, 즉 1089개여야 합니다.
우리는 SparsityPattern 을 셋업해야한다.
임시 구조를 먼저 만들고 0이 아닌 항목에 태그를 지정한 후에
시스템 행렬에서 사용할 수있는 SparsityPattern 개체에 데이터를 복사하여 희소성 패턴을 설정합니다.
DynamicSparsityPattern dsp(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler, dsp);
sparsity_pattern.copy_from(dsp);
SparsityPattern 객체는 행렬의 값을 보유하지 않고 항목이있는 위치 만 저장합니다. 항목 자체는 SparsityPattern 유형의 객체에 저장되며, 그중 system_matrix 변수가 1입니다.
희소성 패턴과 행렬을 구분하여 여러 행렬이 동일한 희소성 패턴을 사용할 수 있도록했습니다. 여기에서는 관련성이 없어 보일 수 있지만 행렬이 가질 수있는 크기와 희소성 패턴을 만드는 데 시간이 걸릴 수 있다는 점을 고려할 때 프로그램에 여러 행렬을 저장해야하는 경우 대규모 문제에서 중요해집니다.
system_matrix.reinit(sparsity_pattern);
이 함수에서 마지막으로 할 일은 우변 벡터와 해 벡터의 크기를 오른쪽 값으로 설정하는 것입니다.
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
}
Step3::assemble_system
다음 단계는 해를 계산하는 선형 시스템을 형성하는 행렬과 우변의 항목을 계산하는 것입니다.
이것은 각 유한 요소 프로그램의 핵심 기능이며 이미 소개의 주요 단계에 대해 논의했습니다.
행렬과 벡터를 어셈블하는 일반적인 접근 방식은 모든 셀을 반복하고 각 셀에서 해당 셀이 전역 행렬에 대한 기여도를 계산하고 구적법(quadrature)으로 오른쪽에 계산하는 것입니다. 지금 깨달아야 할 점은 실제 셀의 구적 점 위치에서 모양 함수의 값이 필요하다는 것입니다. 그러나 유한 요소 모양 함수와 구적 점은 모두 참조 셀에서만 정의됩니다. 따라서 그것들은 우리에게 거의 도움이되지 않으며 실제로 이러한 객체에서 직접 유한 요소 모양 함수 또는 직교 점에 대한 정보를 거의 쿼리하지 않습니다.
오히려 필요한 것은이 데이터를 참조 셀에서 실제 셀로 매핑하는 방법입니다. 이를 수행 할 수있는 클래스는 Mapping에서 파생되지만 다시 한 번 직접 처리 할 필요는 없습니다. 라이브러리의 많은 함수는 매핑 개체를 인수로 사용할 수 있지만 생략하면 단순히 표준 쌍 선형에 의존합니다. Q1 매핑. 우리는이 경로로 이동하고 잠시 신경 쓰지 않을 것입니다 (10 단계, 11 단계 및 12 단계에서 다시 돌아옵니다).
이제 우리가 가지고있는 것은 유한 요소, 구적법, 매핑 객체의 세 가지 클래스 모음입니다. 너무 많기 때문에이 세 가지 간의 정보 교환을 조정하는 클래스 유형이 하나 있습니다. FEValues 클래스입니다. 이러한 개체 3 개 (또는 2 개 및 암시 적 선형 매핑) 각각에 대한 인스턴스가 하나 주어지면 실제 셀의 구적 점에서 모양 함수의 값과 기울기에 대한 정보를 제공 할 수 있습니다.
이 모든 것을 사용하여 다음 함수에서이 문제에 대한 선형 시스템을 조립합니다.
void Step3::assemble_system()
{
QGauss<2> quadrature_formula(fe.degree + 1);
QGauss<2> quadrature_formula(fe.degree + 1);각 셀의 적분을 평가하려면 구적 공식이 필요합니다.
각 방향으로 2 개의 구적 점을 갖는 가우스 공식을 취해 보겠습니다.
이 구적 공식은 최대 3 도의 다항식을 정확히 통합합니다 (1D에서).
이것이 현재 문제에 충분한 지 확인하는 것은 쉽습니다.
FEValues<2> fe_values(fe,
quadrature_formula,
update_values | update_gradients | update_JxW_values);객체를 초기화합니다.
우리가 사용할 유한 요소와 구적 점과 가중치 (Quadrature 객체에 의해 공동으로 설명 됨)를 알아야합니다.
앞서 언급했듯이 우리는 명시 적(explicitly)으로 하나를 지정하는 대신 암시 적(implied) Q1 매핑을 사용합니다.
마지막으로, 우리는 각 셀에서 계산하기를 원하는 것을 말해야합니다.
우리는 구적 점 (우변 (φi, f)의 경우), 기울기 (행렬 항목 ()의 경우)에서 모양 함수의 값이 필요합니다.
(∇φi, ∇φj)), 또한 구적 점의 가중치와 참조 셀에서 실제 셀로의 야 코비 변환의 결정 인자입니다.
우리가 실제로 필요한 정보의이 목록은 FEValues 생성자에 대한 세 번째 인수로 플래그 모음으로 제공됩니다. 이러한 값을 다시 계산하거나 업데이트해야하므로 새 셀로 이동할 때마다 이러한 모든 플래그는 접두사 update_로 시작하고 실제로 업데이트하려는 것이 무엇인지 나타냅니다. 계산 된 모양 함수의 값을 원하는 경우 제공 할 플래그는 update_values;입니다. 그라디언트의 경우 update_gradients입니다. Jacobians의 행렬식과 구적 가중치는 항상 함께 사용되므로 제품 (Jacobians 곱하기 가중치 또는 짧은 JxW) 만 계산됩니다. 필요하기 때문에 update_JxW_values도 나열해야합니다.
이 접근 방식의 장점은 각 셀에 실제로 필요한 정보의 종류를 지정할 수 있다는 것입니다. 이 접근 방식은 필요 여부에 관계없이 각 셀에서 2 차 도함수, 셀에 대한 법선 벡터 등을 포함한 모든 것을 계산하는 접근 방식에 비해 유한 요소 계산 속도를 크게 높일 수 있음을 쉽게 이해할 수 있습니다.
NOTE
구문 update_values | update_gradients | update_JxW_values는
이미 수년 동안 C에서 비트 연산을 프로그래밍하는 데 익숙하지 않은 사람에게 즉시 명확하지 않습니다.
먼저 operator | 는 비트 또는 연산자입니다.
즉, 비트 패턴으로 해석되는 두 개의 정수 인수를 취하고 해당 비트가 두 인수 중 적어도
하나에 설정된 모든 비트가 설정된 정수를 반환합니다. 예를 들어 9 | 10 연산을 생각해보십시오.
이진수에서는 9 = 0b1001 (접두사 0b는 숫자가 이진수로 해석됨을 나타냄) 및 10 = 0b1010입니다.
각 비트를 살펴보고 인수 중 하나에 설정되어 있는지 확인하면
0b1001 | 0b1010 = 0b1011 또는 10 진수 표기법으로 9 | 10 = 11에 도달합니다.
알아야 할 두 번째 정보는 다양한 update_ * 플래그가 정확히 1 비트 세트를 가진 정수라는 것입니다.
예를 들어 update_values = 0b00001 = 1, update_gradients = 0b00010 = 2,
update_JxW_values = 0b10000 = 16이라고 가정합니다.
그런 다음 update_values | update_gradients | update_JxW_values = 0b10011 = 19.
즉, 발생하려는 모든 작업을 나타내는 이진 마스크를 인코딩하는 숫자를 얻습니다.
여기서 각 작업은 정수의 정확히 1 비트에 해당하며 1과 같으면 다음을 의미합니다.
특정 부분은 각 셀에서 업데이트되어야하며, 0이면 계산할 필요가 없음을 의미합니다.
즉, operator | 비트 OR 연산입니다. 이것이 실제로 나타내는 것은
이것과 저것 그리고 다른 것을 원한다는 것입니다.
이러한 바이너리 마스크는 C 프로그래밍에서 매우 일반적이지만
C ++와 같은 상위 수준 언어에서는 그렇지 않을 수 있지만 현재의 목적에 상당히 잘 부합합니다.아래에서 더 사용하기 위해 매우 자주 사용되는 값에 대한 바로 가기를 정의합니다. 즉, 각 셀의 자유도 수에 대한 약어입니다 (2D에 있고 자유도는 꼭지점에만 연결되어 있으므로이 숫자는 4이지만이 변수의 정의를 다음과 같은 방식으로 작성하려고합니다. 나중에 셀당 자유도가 다른 다른 유한 요소를 선택하거나 다른 공간 차원에서 작업하는 것을 배제하지 않습니다.)
일반적으로 이러한 숫자를 알고 있더라도 하드 코딩하는 대신 기호 이름을 사용하는 것이 좋습니다. 예를 들어 유한 요소를 언젠가 변경할 수 있기 때문입니다. 요소 변경은 다른 기능에서 수행되어야하며 프로그램의 다른 부분에서 해당 변경을 수행하는 것을 잊기 쉽습니다. 자신의 계산에 의존하지 않고 대신 올바른 객체에 정보를 요청하는 것이 좋습니다. 여기서는 유한 요소에 셀당 자유도 수에 대해 알려주도록 요청하고 프로그램의 다른 곳에서 선택한 공간 차원 또는 다항식 정도.
여기에있는 단축키는 기본 개념을 설명하기 위해 정의 된 것이지 많은 타이핑을 절약하기 때문이 아니라 다음 루프를 좀 더 읽기 쉽게 만들 것입니다. 큰 프로그램의 여러 위치에서 이러한 바로 가기를 볼 수 있으며 dofs_per_cell은 이러한 종류의 개체에 대한 일반적인 이름입니다.
const unsigned int dofs_per_cell = fe.dofs_per_cell;
이제 우리는 글로벌 행렬과 벡터를 셀 단위로 조합하고 싶다고 말했습니다. 결과를 전역 행렬에 직접 쓸 수는 있지만 희소 행렬의 요소에 대한 액세스가 느리기 때문에 이는 그리 효율적이지 않습니다. 오히려 우리는 먼저 현재 셀에 대한 자유도를 가진 작은 행렬에서 각 셀의 기여도를 계산하고이 셀에 대한 계산이 완료 될 때만 전역 행렬로 전송합니다. 우변 벡터에 대해서도 똑같이합니다. 따라서 먼저 이러한 객체를 할당 해 보겠습니다 (이 객체는 로컬 객체이고 모든 자유도는 다른 모든 자유 도와 결합되며 로컬 연산에 희소 객체가 아닌 전체 행렬 객체를 사용해야합니다. 모든 것은 나중에 전역 희소 행렬로 전송됩니다. 의 위에):
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
각 셀의 기여도를 조합 할 때 자유도의 로컬 번호 (즉, 0에서 dofs_per_cell-1까지의 숫자)를 사용하여이를 수행합니다. 그러나 결과를 전역 행렬로 전송할 때 전역 자유도 수를 알아야합니다. 쿼리 할 때 이러한 숫자에 대한 스크래치 (임시) 배열이 필요합니다 (여기에 사용 된 유형, types :: global_dof_index에 대한 소개 끝에있는 토론 참조).
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
이제 모든 셀에 대한 루프입니다. 우리는 이것이 Triangulation에서 어떻게 작동하는지 전에 보았습니다. DoFHandler 에는 삼각 분할과 정확히 유사한 셀 반복기가 있지만 사용중인 유한 요소의 자유도에 대한 추가 정보가 있습니다. 자유도 핸들러의 활성 셀을 반복하는 것은 Triangulation과 동일하게 작동합니다.
이번에는 셀 유형을 auto 대신 const auto &로 선언합니다. 1 단계에서 세분화 표시기로 플래그를 지정하여 삼각 분할 셀을 수정했습니다. 여기서는 셀을 수정하지 않고 셀만 검사하므로이 불변성을 적용하기 위해 셀을 const로 선언하는 것이 좋습니다.
for (const auto &cell : dof_handler.active_cell_iterators())
{이제 우리는 하나의 셀에 앉아 있으며, 구적 점에서 참조 셀과 실제 셀 간의 매핑에 대한 야 코비 행렬의 결정 요소뿐만 아니라 모양 함수의 값과 기울기를 계산하려고합니다. 이 모든 값은 셀의 지오메트리에 따라 다르기 때문에 FEValues 객체가 각 셀에서 다시 계산해야합니다.
fe_values.reinit(cell);Next, reset the local cell's contributions to global matrix and global right hand side to zero, before we fill them:
cell_matrix = 0;
cell_rhs = 0;
이제 셀에 대한 통합을 시작할 때입니다. 모든 구적 점을 반복하여 수행합니다.이 작업은 q_index로 번호를 매 깁니다.
for (const unsigned int q_index : fe_values.quadrature_point_indices())
{먼저 행렬을 조립합니다. 라플라스 문제의 경우 각 셀의 행렬은 모양 함수 i 및 j의 기울기에 대한 적분입니다. 적분하지 않고 구적 점을 사용하기 때문에 이것은 적분의 모든 구적 점에 대한 합과 구적 점에서 야 코비 행렬의 행렬식 곱하기이 구적 점의 가중치를 곱한 값입니다. fe_values.shape_grad (i, q_index);를 사용하여 숫자 q_index를 가진 구적 점에서 모양 함수 i의 기울기를 얻을 수 있습니다. 이 그래디언트는 2 차원 벡터 (실제로는 Tensor <1, dim> 유형이며 여기서는 dim = 2 임)이며 이러한 두 벡터의 곱은 스칼라 곱입니다. 즉, 두 shape_grad 함수 호출의 곱은 다음과 같습니다. 내적. 이것은 차례로 Jacobian determinant와 quadrature point weight (FEValues::JxW()에 대한 호출로 합쳐 짐)를 곱합니다. 마지막으로, 이것은 모든 모양 함수 i와 j에 대해 반복됩니다.
for (const unsigned int i : fe_values.dof_indices())
for (const unsigned int j : fe_values.dof_indices())
cell_matrix(i, j) +=
(fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
fe_values.JxW(q_index)); // dx
그런 다음 오른쪽에 대해서도 동일한 작업을 수행합니다. 여기서 적분은 형상 함수 i 곱하기 우변 함수이며 상수 값이 1 인 함수로 선택합니다 (다음 프로그램에서 더 흥미로운 예를 고려할 것입니다).
for (const unsigned int i : fe_values.dof_indices())
cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
1 * // f(x_q)
fe_values.JxW(q_index)); // dx
}
이제이 셀의 기여도를 얻었으므로 이를 전역 행렬과 오른쪽으로 전송해야합니다.
이를 위해 우리는 먼저이 세포의 자유도가 어떤 글로벌 숫자인지 알아 내야합니다.
셀에 해당 정보를 요청 해 보겠습니다.
cell->get_dof_indices(local_dof_indices);
그런 다음 모든 모양 함수 i 및 j를 다시 반복하고 로컬 요소를 전역 행렬로 전송합니다. 전역 번호는 local_dof_indices [i]를 사용하여 얻을 수 있습니다.
for (const unsigned int i : fe_values.dof_indices())
for (const unsigned int j : fe_values.dof_indices())
system_matrix.add(local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i, j));
And again, we do the same thing for the right hand side vector.
for (const unsigned int i : fe_values.dof_indices())
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
이제 거의 모든 것이 개별 시스템의 솔루션을 위해 설정되었습니다. 그러나 우리는 아직 경계 값을 처리하지 않았습니다 (사실 Dirichlet 경계 값이없는 라플라스 방정식은 이산 솔루션에 임의의 상수를 추가 할 수 있기 때문에 고유하게 풀 수 없습니다). 그러므로 우리는 상황에 대해 뭔가를해야합니다.
이를 위해 먼저 경계의 자유도 목록과 형상 함수가 가질 값을 얻습니다. 단순화를 위해 경계 값 함수를 경계에 투영하는 대신 보간 만합니다. 라이브러리에 정확히 다음과 같은 함수가 있습니다 : VectorTools::interpolate_boundary_values(). 매개 변수는 다음과 같습니다 (기본값이 존재하고 우리가 신경 쓰지 않는 매개 변수 생략). 경계에서 자유도의 전역 수를 가져 오는 DoFHandler 객체; 경계 값이 보간되는 경계의 구성 요소 경계 값 기능 자체; 및 출력 개체.
경계의 구성 요소는 다음과 같습니다. 많은 경우 경계의 일부에만 특정 경계 값을 적용 할 수 있습니다. 예를 들어 유체 역학에서 유입 및 유출 경계가 있거나 바디의 변형 계산에서 바디의 고정 및 자유 부분이있을 수 있습니다. 그런 다음 경계의 서로 다른 부분을 표시기로 표시하고 interpolate_boundary_values 함수에 경계의 특정 부분 (예 : 고정 된 부분 또는 유입 경계)에서만 경계 값을 계산하도록 지시 할 수 있습니다. 기본적으로 모든 경계에는 달리 지정하지 않는 한 0 경계 표시기가 있습니다. 경계 섹션에 다른 경계 조건이있는 경우 다른 경계 표시기로 해당 부품에 번호를 지정해야합니다. 그러면 아래의 함수 호출은 경계 표시기가 실제로 두 번째 인수로 지정된 0 인 경계 부분에 대한 경계 값만 결정합니다.
경계 값을 설명하는 함수는 Function 유형 또는 파생 클래스의 개체입니다. 파생 클래스 중 하나는 Functions::ZeroFunction으로, 모든 곳에서 0 인 함수를 (예기치 않게) 설명합니다. 그런 객체를 제자리에 만들어 VectorTools::interpolate_boundary_values() 함수에 전달합니다.
마지막으로 출력 객체는 전역 자유도 수 (즉, 경계의 자유도 수)와 해당 경계 값 (여기서는 모든 항목에 대해 0 임) 쌍의 목록입니다. DoF 번호를 경계 값에 매핑하는 것은 std :: map 클래스에 의해 수행됩니다.
std::map<types::global_dof_index, double> boundary_values;
VectorTools::interpolate_boundary_values(dof_handler,
0,
Functions::ZeroFunction<2>(),
boundary_values);이제 경계 DoF 목록과 각각의 경계 값을 얻었으므로이를 사용하여
방정식 시스템을 적절하게 수정 해 보겠습니다. 이것은 다음 함수 호출에 의해 수행됩니다.
MatrixTools::apply_boundary_values(boundary_values,
system_matrix,
solution,
system_rhs);
}
Step3::solve
다음 함수는 이산화 방정식을 간단히 해결합니다. 가우스 제거 또는 LU 분해와 같은 직접 솔버의 경우 시스템이 상당히 크기 때문에 Conjugate Gradient 알고리즘을 사용합니다. 여기서 변수의 수 (단지 1089)는 유한 요소 계산의 경우 매우 적은 수이며 여기서 100.000은 더 일반적인 수입니다. 이 변수 수에 대해 직접 메서드를 더 이상 사용할 수 없으며 CG와 같은 메서드를 사용해야합니다.
void Step3::solve()
{
먼저, CG 알고리즘에 언제 멈출 지 알려주는 방법을 아는 객체가 필요합니다. 이 작업은 SolverControl 객체를 사용하여 수행되며 중지 기준으로 다음과 같이 말합니다. 최대 1000 회 반복 후 중지 (1089 개의 변수에 필요한 것보다 훨씬 더 많습니다. 실제로 사용 된 변수는 결과 섹션 참조) 잔차의 노름이 10-12 미만이면 중지합니다. 실제로 후자의 기준은 반복을 중지하는 기준이됩니다.
SolverControl solver_control(1000, 1e-12);
그런 다음 우리는 솔버 자체가 필요합니다. SolverCG 클래스의 템플릿 매개 변수는 벡터의 유형이지만 빈 꺾쇠 괄호는 단순히 기본 인수 (Vector <double>)를 사용함을 나타냅니다.
SolverCG<Vector<double>> solver(solver_control);
이제 연립 방정식을 풉니 다. CG 솔버는 선 조건자를 네 번째 인수로 사용합니다. 우리는 아직 이것에 대해 탐구 할 준비가되어 있지 않다고 느낍니다. 그래서 우리는 신원 연산을 선 조건 자로 사용하라고 말합니다.
solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());이제 솔버가 작업을 완료 했으므로 솔루션 변수에는 솔루션 함수의 절점 값이 포함됩니다.
}
Step3::output_results
일반적인 유한 요소 프로그램의 마지막 부분은 결과를 출력하고 일부 후 처리를 수행하는 것입니다
(예 : 경계에서 최대 응력 값 또는 유출에 걸친 평균 플럭스 계산 등).
여기에는 그러한 후 처리가 없지만 솔루션을 파일에 작성하고 싶습니다.
void Step3::output_results() const
{
출력을 파일에 쓰려면 출력 형식 등에 대해 알고있는 객체가 필요합니다. 이것이 DataOut 클래스이며 해당 유형의 객체가 필요합니다.
DataOut<2> data_out;
이제 우리는 그것이 쓸 값을 어디에서 가져올 것인지 알려야합니다. 사용할 DoFHandler 객체와 솔루션 벡터 (및 솔루션 변수가 출력 파일에 나타날 이름)를 알려줍니다. 출력에서보고 싶은 벡터가 두 개 이상있는 경우 (예 : 오른쪽, 셀당 오류 등) 다음과 같이 추가합니다.
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "solution");DataOut 객체가 작업 할 데이터를 알고 나면 백엔드가 처리 할 수있는 것으로 처리하도록 지정해야합니다. 그 이유는 프런트 엔드 (DoFHandler 객체 및 데이터 벡터를 처리하는 방법에 대해 알고 있음)를 백엔드 (많은 다른 출력 형식을 알고 있음)에서 분리하고 중간 데이터 형식을 사용하여 프런트 엔드에서 백엔드로 데이터를 전송하기 때문입니다. . 데이터는 다음 함수에 의해이 중간 형식으로 변환됩니다.
data_out.build_patches();이제 실제 출력을위한 모든 것이 준비되었습니다. VTK 형식을 사용하여 파일을 열고 여기에 데이터를 쓰십시오 (여기에서 사용중인 DataOut 클래스에는 포스트 스크립트, AVS, GMV, Gnuplot 또는 기타 파일 형식으로 데이터를 쓸 수있는 다른 많은 함수가 있습니다) :
std::ofstream output("solution.vtk");
data_out.write_vtk(output);
}
Step3::run
마지막으로이 클래스의 마지막 함수는 Step3 클래스의 다른 모든 함수를 호출하는 주 함수입니다. 이것이 수행되는 순서는 대부분의 유한 요소 프로그램이 작동하는 순서와 유사합니다. 이름은 대부분 자명하기 때문에 설명 할 내용이 많지 않습니다.
void Step3::run()
{
make_grid();
setup_system();
assemble_system();
solve();
output_results();
}
Main function
이것이 프로그램의 main function입니다. 메인 함수의 개념은 대부분 C ++ 프로그래밍 이전의 객체 지향 이전 시대의 잔재이기 때문에 종종 최상위 클래스의 객체를 생성하고 기본 함수를 호출하는 것 이상을 수행하지 않습니다.
마지막으로, 함수의 첫 번째 줄은 deal.II가 생성 할 수있는 일부 진단의 출력을 활성화하는 데 사용됩니다. deallog 변수 (deallog가 아니라 deal-log를 나타냄)는 라이브러리의 일부가 출력을 기록하는 스트림을 나타냅니다. 예를 들어, 반복 솔버는이 튜토리얼 프로그램을 실행할 때 볼 수있는 진단 (시작 잔차, 솔버 단계 수, 최종 잔차)을 생성합니다.
deallog의 출력은 콘솔, 파일 또는 둘 다에 기록 될 수 있습니다. 수년에 걸쳐 사용자가 명시 적으로 요청할 때만 프로그램이 출력을 생성해야한다는 것을 배웠기 때문에 둘 다 기본적으로 비활성화되어 있습니다. 그러나 이것은 변경 될 수 있으며 이것이 수행되는 방법을 설명하기 위해 deallog의 작동 방식을 설명해야합니다. 라이브러리의 개별 부분이 출력을 기록하려고 할 때이 출력이 들어갈 "컨텍스트"또는 "섹션"을 엽니 다. 배치. 출력을 작성하려는 부분의 끝에서이 섹션을 다시 종료합니다. 함수는이 출력 섹션이 열려있는 범위 내에서 다른 함수를 호출 할 수 있으므로 실제로 출력은 이러한 섹션에 계층 적으로 중첩 될 수 있습니다. deallog가 변수 인 LogStream 클래스는 모든 출력이 콜론으로 구분 된 접두사와 함께 행의 왼쪽 끝에이 접두사와 함께 인쇄되기 때문에 이러한 각 섹션을 "접두사"라고 부릅니다. 항상 "DEAL"이라는 기본 접두사가 있습니다 ( "DEAL"이라는 이전 라이브러리의 후속 버전 인 deal.II의 역사에 대한 힌트이며 여기서 LogStream 클래스는 거래에 사용 된 몇 안되는 코드 중 하나입니다 .II). ).
기본적으로 logstream은 접두사가 0 인 행만 출력합니다. 즉, 기본 "DEAL"접두사가 항상 있기 때문에 모든 출력이 비활성화됩니다. 그러나 출력되어야하는 라인에 대해 다른 최대 접두사 수를 설정할 수 있으며 실제로 여기서는 LogStream::depth_console()을 호출하여 두 개로 설정합니다. 즉, 모든 화면 출력에 대해 기본 "DEAL"을 넘어서 하나의 추가 접두사를 푸시 한 컨텍스트는 해당 출력을 화면 ( "콘솔")에 인쇄 할 수있는 반면, 세 개 이상의 접두사가 활성화 된 모든 추가 중첩 섹션 deallog에 기록하지만 deallog는이 출력을 화면에 전달하지 않습니다. 따라서이 예제를 실행하거나 "결과"섹션을 보면 "DEAL : CG"접두사가 붙은 솔버 통계가 표시됩니다. 이것은 현재 프로그램의 컨텍스트에 충분하지만 나중에 솔버가 더 깊게 중첩되고 깊이를 더 높게 설정하여 유용한 정보를 얻을 수있는 예를 나중에 볼 수 있습니다 (예 : 22 단계).
int main()
{
deallog.depth_console(2);
Step3 laplace_problem;
laplace_problem.run();
return 0;
}
Results
The output of the program looks as follows:
Number of active cells: 1024
Number of degrees of freedom: 1089
DEAL:cg::Starting value 0.121094
DEAL:cg::Convergence step 48 value 5.33692e-13The first two lines is what we wrote to cout. The last two lines were generated without our intervention by the CG solver. The first two lines state the residual at the start of the iteration, while the last line tells us that the solver needed 47 iterations to bring the norm of the residual to 5.3e-13, i.e. below the threshold 1e-12 which we have set in the `solve' function. We will show in the next program how to suppress this output, which is sometimes useful for debugging purposes, but often clutters up the screen display.
Apart from the output shown above, the program generated the file solution.vtk, which is in the VTK format that is widely used by many visualization programs today – including the two heavy-weights VisIt and Paraview that are the most commonly used programs for this purpose today.
Using VisIt, it is not very difficult to generate a picture of the solution like this:
#include <deal.II/grid/tria.h>
#include <deal.II/dofs/dof_handler.h>
// And this is the file in which the functions are declared that create grids:
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.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>
// Finally, this is for output to a file and to the console:
#include <deal.II/numerics/data_out.h>
#include <fstream>
#include <iostream>
// ...and this is to import the deal.II namespace into the global scope:
using namespace dealii;
// @sect3{The <code>Step3</code> class}
class Step3
{
public:
Step3();
void run();
private:
void make_grid();
void setup_system();
void assemble_system();
void solve();
void output_results() const;
Triangulation<3> triangulation;
FE_Q<3> fe;
DoFHandler<3> dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};
// @sect4{Step3::Step3}
Step3::Step3()
: fe(1)
, dof_handler(triangulation)
{}
// @sect4{Step3::make_grid}
void Step3::make_grid()
{
GridGenerator::hyper_cube(triangulation, -1, 1);
triangulation.refine_global(5);
std::cout << "Number of active cells: " << triangulation.n_active_cells()
<< std::endl;
}
// @sect4{Step3::setup_system}
void Step3::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());
}
// @sect4{Step3::assemble_system}
void Step3::assemble_system()
{
QGauss<3> quadrature_formula(fe.degree + 1);
FEValues<3> fe_values(fe,
quadrature_formula,
update_values | update_gradients | 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);
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
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)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
cell_matrix(i, j) +=
(fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
fe_values.JxW(q_index)); // dx
for (unsigned int i = 0; i < dofs_per_cell; ++i)
cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
1 * // f(x_q)
fe_values.JxW(q_index)); // dx
}
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)
system_matrix.add(local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i, j));
// And again, we do the same thing for the right hand side vector.
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,
Functions::ZeroFunction<3>(),
boundary_values);
MatrixTools::apply_boundary_values(boundary_values,
system_matrix,
solution,
system_rhs);
}
// @sect4{Step3::solve}
void Step3::solve()
{
SolverControl solver_control(1000, 1e-12);
SolverCG<> solver(solver_control);
solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
}
// @sect4{Step3::output_results}
void Step3::output_results() const
{
DataOut<3> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "solution");
data_out.build_patches();
std::ofstream output("solution.vtk");
data_out.write_vtk(output);
}
// @sect4{Step3::run}
void Step3::run()
{
make_grid();
setup_system();
assemble_system();
solve();
output_results();
}
// @sect3{The <code>main</code> function}
int main()
{
deallog.depth_console(2);
Step3 laplace_problem;
laplace_problem.run();
return 0;
}
/* ---------------------------------------------------------------------
*
* Copyright (C) 1999 - 2019 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.md at
* the top level directory of deal.II.
*
* ---------------------------------------------------------------------
*
* Authors: Wolfgang Bangerth, 1999,
* Guido Kanschat, 2011
*/
// These include files are already known to you. They declare the classes
// which handle triangulations and enumeration of degrees of freedom:
#include <deal.II/grid/tria.h>
#include <deal.II/dofs/dof_handler.h>
// And this is the file in which the functions are declared that create grids:
#include <deal.II/grid/grid_generator.h>
// The next three files contain classes which are needed for loops over all
// cells and to get the information from the cell objects. The first two have
// been used before to get geometric information from cells; the last one is
// new and provides information about the degrees of freedom local to a cell:
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_accessor.h>
// This file contains the description of the Lagrange interpolation finite
// element:
#include <deal.II/fe/fe_q.h>
// And this file is needed for the creation of sparsity patterns of sparse
// matrices, as shown in previous examples:
#include <deal.II/dofs/dof_tools.h>
// The next two files are needed for assembling the matrix using quadrature on
// each cell. The classes declared in them will be explained below:
#include <deal.II/fe/fe_values.h>
#include <deal.II/base/quadrature_lib.h>
// The following three include files we need for the treatment of boundary
// values:
#include <deal.II/base/function.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
// We're now almost to the end. The second to last group of include files is
// for the linear algebra which we employ to solve the system of equations
// arising from the finite element discretization of the Laplace equation. We
// will use vectors and full matrices for assembling the system of equations
// locally on each cell, and transfer the results into a sparse matrix. We
// will then use a Conjugate Gradient solver to solve the problem, for which
// we need a preconditioner (in this program, we use the identity
// preconditioner which does nothing, but we need to include the file anyway):
#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>
// Finally, this is for output to a file and to the console:
#include <deal.II/numerics/data_out.h>
#include <fstream>
#include <iostream>
// ...and this is to import the deal.II namespace into the global scope:
using namespace dealii;
// @sect3{The <code>Step3</code> class}
// Instead of the procedural programming of previous examples, we encapsulate
// everything into a class for this program. The class consists of functions
// which each perform certain aspects of a finite element program, a `main'
// function which controls what is done first and what is done next, and a
// list of member variables.
// The public part of the class is rather short: it has a constructor and a
// function `run' that is called from the outside and acts as something like
// the `main' function: it coordinates which operations of this class shall be
// run in which order. Everything else in the class, i.e. all the functions
// that actually do anything, are in the private section of the class:
class Step3
{
public:
Step3();
void run();
// Then there are the member functions that mostly do what their names
// suggest and whose have been discussed in the introduction already. Since
// they do not need to be called from outside, they are made private to this
// class.
private:
void make_grid();
void setup_system();
void assemble_system();
void solve();
void output_results() const;
// And finally we have some member variables. There are variables describing
// the triangulation and the global numbering of the degrees of freedom (we
// will specify the exact polynomial degree of the finite element in the
// constructor of this class)...
Triangulation<3> triangulation;
FE_Q<3> fe;
DoFHandler<3> dof_handler;
// ...variables for the sparsity pattern and values of the system matrix
// resulting from the discretization of the Laplace equation...
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
// ...and variables which will hold the right hand side and solution
// vectors.
Vector<double> solution;
Vector<double> system_rhs;
};
// @sect4{Step3::Step3}
// Here comes the constructor. It does not much more than first to specify
// that we want bi-linear elements (denoted by the parameter to the finite
// element object, which indicates the polynomial degree), and to associate
// the dof_handler variable to the triangulation we use. (Note that the
// triangulation isn't set up with a mesh at all at the present time, but the
// DoFHandler doesn't care: it only wants to know which triangulation it will
// be associated with, and it only starts to care about an actual mesh once
// you try to distribute degree of freedom on the mesh using the
// distribute_dofs() function.) All the other member variables of the Step3
// class have a default constructor which does all we want.
Step3::Step3()
: fe(1)
, dof_handler(triangulation)
{}
// @sect4{Step3::make_grid}
// Now, the first thing we've got to do is to generate the triangulation on
// which we would like to do our computation and number each vertex with a
// degree of freedom. We have seen these two steps in step-1 and step-2
// before, respectively.
//
// This function does the first part, creating the mesh. We create the grid
// and refine all cells five times. Since the initial grid (which is the
// square $[-1,1] \times [-1,1]$) consists of only one cell, the final grid
// has 32 times 32 cells, for a total of 1024.
//
// Unsure that 1024 is the correct number? We can check that by outputting the
// number of cells using the <code>n_active_cells()</code> function on the
// triangulation.
void Step3::make_grid()
{
GridGenerator::hyper_cube(triangulation, -1, 1);
triangulation.refine_global(5);
std::cout << "Number of active cells: " << triangulation.n_active_cells()
<< std::endl;
}
// @note We call the Triangulation::n_active_cells() function, rather than
// Triangulation::n_cells(). Here, <i>active</i> means the cells that aren't
// refined any further. We stress the adjective "active" since there are more
// cells, namely the parent cells of the finest cells, their parents, etc, up
// to the one cell which made up the initial grid. Of course, on the next
// coarser level, the number of cells is one quarter that of the cells on the
// finest level, i.e. 256, then 64, 16, 4, and 1. If you called
// <code>triangulation.n_cells()</code> instead in the code above, you would
// consequently get a value of 1365 instead. On the other hand, the number of
// cells (as opposed to the number of active cells) is not typically of much
// interest, so there is no good reason to print it.
// @sect4{Step3::setup_system}
// Next we enumerate all the degrees of freedom and set up matrix and vector
// objects to hold the system data. Enumerating is done by using
// DoFHandler::distribute_dofs(), as we have seen in the step-2 example. Since
// we use the FE_Q class and have set the polynomial degree to 1 in the
// constructor, i.e. bilinear elements, this associates one degree of freedom
// with each vertex. While we're at generating output, let us also take a look
// at how many degrees of freedom are generated:
void Step3::setup_system()
{
dof_handler.distribute_dofs(fe);
std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
// There should be one DoF for each vertex. Since we have a 32 times 32
// grid, the number of DoFs should be 33 times 33, or 1089.
// As we have seen in the previous example, we set up a sparsity pattern by
// first creating a temporary structure, tagging those entries that might be
// nonzero, and then copying the data over to the SparsityPattern object
// that can then be used by the system matrix.
DynamicSparsityPattern dsp(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler, dsp);
sparsity_pattern.copy_from(dsp);
// Note that the SparsityPattern object does not hold the values of the
// matrix, it only stores the places where entries are. The entries
// themselves are stored in objects of type SparseMatrix, of which our
// variable system_matrix is one.
//
// The distinction between sparsity pattern and matrix was made to allow
// several matrices to use the same sparsity pattern. This may not seem
// relevant here, but when you consider the size which matrices can have,
// and that it may take some time to build the sparsity pattern, this
// becomes important in large-scale problems if you have to store several
// matrices in your program.
system_matrix.reinit(sparsity_pattern);
// The last thing to do in this function is to set the sizes of the right
// hand side vector and the solution vector to the right values:
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
}
// @sect4{Step3::assemble_system}
// The next step is to compute the entries of the matrix and right hand side
// that form the linear system from which we compute the solution. This is the
// central function of each finite element program and we have discussed the
// primary steps in the introduction already.
//
// The general approach to assemble matrices and vectors is to loop over all
// cells, and on each cell compute the contribution of that cell to the global
// matrix and right hand side by quadrature. The point to realize now is that
// we need the values of the shape functions at the locations of quadrature
// points on the real cell. However, both the finite element shape functions
// as well as the quadrature points are only defined on the reference
// cell. They are therefore of little help to us, and we will in fact hardly
// ever query information about finite element shape functions or quadrature
// points from these objects directly.
//
// Rather, what is required is a way to map this data from the reference cell
// to the real cell. Classes that can do that are derived from the Mapping
// class, though one again often does not have to deal with them directly:
// many functions in the library can take a mapping object as argument, but
// when it is omitted they simply resort to the standard bilinear Q1
// mapping. We will go this route, and not bother with it for the moment (we
// come back to this in step-10, step-11, and step-12).
//
// So what we now have is a collection of three classes to deal with: finite
// element, quadrature, and mapping objects. That's too much, so there is one
// type of class that orchestrates information exchange between these three:
// the FEValues class. If given one instance of each three of these objects
// (or two, and an implicit linear mapping), it will be able to provide you
// with information about values and gradients of shape functions at
// quadrature points on a real cell.
//
// Using all this, we will assemble the linear system for this problem in the
// following function:
void Step3::assemble_system()
{
// Ok, let's start: we need a quadrature formula for the evaluation of the
// integrals on each cell. Let's take a Gauss formula with two quadrature
// points in each direction, i.e. a total of four points since we are in
// 2D. This quadrature formula integrates polynomials of degrees up to three
// exactly (in 1D). It is easy to check that this is sufficient for the
// present problem:
QGauss<3> quadrature_formula(fe.degree + 1);
// And we initialize the object which we have briefly talked about above. It
// needs to be told which finite element we want to use, and the quadrature
// points and their weights (jointly described by a Quadrature object). As
// mentioned, we use the implied Q1 mapping, rather than specifying one
// ourselves explicitly. Finally, we have to tell it what we want it to
// compute on each cell: we need the values of the shape functions at the
// quadrature points (for the right hand side $(\varphi_i,f)$), their
// gradients (for the matrix entries $(\nabla \varphi_i, \nabla
// \varphi_j)$), and also the weights of the quadrature points and the
// determinants of the Jacobian transformations from the reference cell to
// the real cells.
//
// This list of what kind of information we actually need is given as a
// collection of flags as the third argument to the constructor of
// FEValues. Since these values have to be recomputed, or updated, every
// time we go to a new cell, all of these flags start with the prefix
// <code>update_</code> and then indicate what it actually is that we want
// updated. The flag to give if we want the values of the shape functions
// computed is #update_values; for the gradients it is
// #update_gradients. The determinants of the Jacobians and the quadrature
// weights are always used together, so only the products (Jacobians times
// weights, or short <code>JxW</code>) are computed; since we need them, we
// have to list #update_JxW_values as well:
FEValues<3> fe_values(fe,
quadrature_formula,
update_values | update_gradients | update_JxW_values);
// The advantage of this approach is that we can specify what kind of
// information we actually need on each cell. It is easily understandable
// that this approach can significantly speed up finite element computations,
// compared to approaches where everything, including second derivatives,
// normal vectors to cells, etc are computed on each cell, regardless of
// whether they are needed or not.
//
// @note The syntax <code>update_values | update_gradients |
// update_JxW_values</code> is not immediately obvious to anyone not
// used to programming bit operations in C for years already. First,
// <code>operator|</code> is the <i>bitwise or operator</i>, i.e.,
// it takes two integer arguments that are interpreted as bit
// patterns and returns an integer in which every bit is set for
// which the corresponding bit is set in at least one of the two
// arguments. For example, consider the operation
// <code>9|10</code>. In binary, <code>9=0b1001</code> (where the
// prefix <code>0b</code> indicates that the number is to be
// interpreted as a binary number) and <code>10=0b1010</code>. Going
// through each bit and seeing whether it is set in one of the
// argument, we arrive at <code>0b1001|0b1010=0b1011</code> or, in
// decimal notation, <code>9|10=11</code>. The second piece of
// information you need to know is that the various
// <code>update_*</code> flags are all integers that have <i>exactly
// one bit set</i>. For example, assume that
// <code>update_values=0b00001=1</code>,
// <code>update_gradients=0b00010=2</code>,
// <code>update_JxW_values=0b10000=16</code>. Then
// <code>update_values | update_gradients | update_JxW_values =
// 0b10011 = 19</code>. In other words, we obtain a number that
// <i>encodes a binary mask representing all of the operations you
// want to happen</i>, where each operation corresponds to exactly
// one bit in the integer that, if equal to one, means that a
// particular piece should be updated on each cell and, if it is
// zero, means that we need not compute it. In other words, even
// though <code>operator|</code> is the <i>bitwise OR operation</i>,
// what it really represents is <i>I want this AND that AND the
// other</i>. Such binary masks are quite common in C programming,
// but maybe not so in higher level languages like C++, but serve
// the current purpose quite well.
// For use further down below, we define two shortcuts for values that will
// be used very frequently. First, an abbreviation for the number of degrees
// of freedom on each cell (since we are in 2D and degrees of freedom are
// associated with vertices only, this number is four, but we rather want to
// write the definition of this variable in a way that does not preclude us
// from later choosing a different finite element that has a different
// number of degrees of freedom per cell, or work in a different space
// dimension).
//
// Secondly, we also define an abbreviation for the number of quadrature
// points (here that should be four). In general, it is a good idea to use
// their symbolic names instead of hard-coding these numbers even if you know
// them, since you may want to change the quadrature formula and/or finite
// element at some time; the program will just work with these changes,
// without the need to change anything in this function.
//
// The shortcuts, finally, are only defined to make the following loops a
// bit more readable. You will see them in many places in larger programs,
// and `dofs_per_cell` and `n_q_points` are more or less by convention the
// standard names for these purposes:
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
// Now, we said that we wanted to assemble the global matrix and vector
// cell-by-cell. We could write the results directly into the global matrix,
// but this is not very efficient since access to the elements of a sparse
// matrix is slow. Rather, we first compute the contribution of each cell in
// a small matrix with the degrees of freedom on the present cell, and only
// transfer them to the global matrix when the computations are finished for
// this cell. We do the same for the right hand side vector. So let's first
// allocate these objects (these being local objects, all degrees of freedom
// are coupling with all others, and we should use a full matrix object
// rather than a sparse one for the local operations; everything will be
// transferred to a global sparse matrix later on):
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
// When assembling the contributions of each cell, we do this with the local
// numbering of the degrees of freedom (i.e. the number running from zero
// through dofs_per_cell-1). However, when we transfer the result into the
// global matrix, we have to know the global numbers of the degrees of
// freedom. When we query them, we need a scratch (temporary) array for
// these numbers (see the discussion at the end of the introduction for
// the type, types::global_dof_index, used here):
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
// Now for the loop over all cells. We have seen before how this works for a
// triangulation. A DoFHandler has cell iterators that are exactly analogous
// to those of a Triangulation, but with extra information about the degrees
// of freedom for the finite element you're using. Looping over the active
// cells of a degree-of-freedom handler works the same as for a triangulation.
//
// Note that we declare the type of the cell as `const auto &` instead of
// `auto` this time around. In step 1, we were modifying the cells of the
// triangulation by flagging them with refinement indicators. Here we're only
// examining the cells without modifying them, so it's good practice to
// declare `cell` as `const` in order to enforce this invariant.
for (const auto &cell : dof_handler.active_cell_iterators())
{
// We are now sitting on one cell, and we would like the values and
// gradients of the shape functions be computed, as well as the
// determinants of the Jacobian matrices of the mapping between
// reference cell and true cell, at the quadrature points. Since all
// these values depend on the geometry of the cell, we have to have the
// FEValues object re-compute them on each cell:
fe_values.reinit(cell);
// Next, reset the local cell's contributions to global matrix and
// global right hand side to zero, before we fill them:
cell_matrix = 0;
cell_rhs = 0;
// Now it is time to start integration over the cell, which we
// do by looping over all quadrature points, which we will
// number by q_index.
for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
{
// First assemble the matrix: For the Laplace problem, the
// matrix on each cell is the integral over the gradients of
// shape function i and j. Since we do not integrate, but
// rather use quadrature, this is the sum over all
// quadrature points of the integrands times the determinant
// of the Jacobian matrix at the quadrature point times the
// weight of this quadrature point. You can get the gradient
// of shape function $i$ at quadrature point with number q_index by
// using <code>fe_values.shape_grad(i,q_index)</code>; this
// gradient is a 2-dimensional vector (in fact it is of type
// Tensor@<1,dim@>, with here dim=2) and the product of two
// such vectors is the scalar product, i.e. the product of
// the two shape_grad function calls is the dot
// product. This is in turn multiplied by the Jacobian
// determinant and the quadrature point weight (that one
// gets together by the call to FEValues::JxW() ). Finally,
// this is repeated for all shape functions $i$ and $j$:
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) * // grad phi_i(x_q)
fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
fe_values.JxW(q_index)); // dx
// We then do the same thing for the right hand side. Here,
// the integral is over the shape function i times the right
// hand side function, which we choose to be the function
// with constant value one (more interesting examples will
// be considered in the following programs).
for (unsigned int i = 0; i < dofs_per_cell; ++i)
cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
1 * // f(x_q)
fe_values.JxW(q_index)); // dx
}
// Now that we have the contribution of this cell, we have to transfer
// it to the global matrix and right hand side. To this end, we first
// have to find out which global numbers the degrees of freedom on this
// cell have. Let's simply ask the cell for that information:
cell->get_dof_indices(local_dof_indices);
// Then again loop over all shape functions i and j and transfer the
// local elements to the global matrix. The global numbers can be
// obtained using local_dof_indices[i]:
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
system_matrix.add(local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i, j));
// And again, we do the same thing for the right hand side vector.
for (unsigned int i = 0; i < dofs_per_cell; ++i)
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
// Now almost everything is set up for the solution of the discrete
// system. However, we have not yet taken care of boundary values (in fact,
// Laplace's equation without Dirichlet boundary values is not even uniquely
// solvable, since you can add an arbitrary constant to the discrete
// solution). We therefore have to do something about the situation.
//
// For this, we first obtain a list of the degrees of freedom on the
// boundary and the value the shape function shall have there. For
// simplicity, we only interpolate the boundary value function, rather than
// projecting it onto the boundary. There is a function in the library which
// does exactly this: VectorTools::interpolate_boundary_values(). Its
// parameters are (omitting parameters for which default values exist and
// that we don't care about): the DoFHandler object to get the global
// numbers of the degrees of freedom on the boundary; the component of the
// boundary where the boundary values shall be interpolated; the boundary
// value function itself; and the output object.
//
// The component of the boundary is meant as follows: in many cases, you may
// want to impose certain boundary values only on parts of the boundary. For
// example, you may have inflow and outflow boundaries in fluid dynamics, or
// clamped and free parts of bodies in deformation computations of
// bodies. Then you will want to denote these different parts of the
// boundary by different numbers and tell the interpolate_boundary_values
// function to only compute the boundary values on a certain part of the
// boundary (e.g. the clamped part, or the inflow boundary). By default, all
// boundaries have the number `0`, and since we have not changed that, this
// is still so; therefore, if we give `0` as the desired portion of the
// boundary, this means we get the whole boundary. If you have boundaries
// with kinds of boundaries, you have to number them differently. The
// function call below will then only determine boundary values for parts of
// the boundary.
//
// The function describing the boundary values is an object of type Function
// or of a derived class. One of the derived classes is
// Functions::ZeroFunction, which describes (not unexpectedly) a function
// which is zero everywhere. We create such an object in-place and pass it to
// the VectorTools::interpolate_boundary_values() function.
//
// Finally, the output object is a list of pairs of global degree of freedom
// numbers (i.e. the number of the degrees of freedom on the boundary) and
// their boundary values (which are zero here for all entries). This mapping
// of DoF numbers to boundary values is done by the <code>std::map</code>
// class.
std::map<types::global_dof_index, double> boundary_values;
VectorTools::interpolate_boundary_values(dof_handler,
0,
Functions::ZeroFunction<3>(),
boundary_values);
// Now that we got the list of boundary DoFs and their respective boundary
// values, let's use them to modify the system of equations
// accordingly. This is done by the following function call:
MatrixTools::apply_boundary_values(boundary_values,
system_matrix,
solution,
system_rhs);
}
// @sect4{Step3::solve}
// The following function simply solves the discretized equation. As the
// system is quite a large one for direct solvers such as Gauss elimination or
// LU decomposition, we use a Conjugate Gradient algorithm. You should
// remember that the number of variables here (only 1089) is a very small
// number for finite element computations, where 100.000 is a more usual
// number. For this number of variables, direct methods are no longer usable
// and you are forced to use methods like CG.
void Step3::solve()
{
// First, we need to have an object that knows how to tell the CG algorithm
// when to stop. This is done by using a SolverControl object, and as
// stopping criterion we say: stop after a maximum of 1000 iterations (which
// is far more than is needed for 1089 variables; see the results section to
// find out how many were really used), and stop if the norm of the residual
// is below $10^{-12}$. In practice, the latter criterion will be the one
// which stops the iteration:
SolverControl solver_control(1000, 1e-12);
// Then we need the solver itself. The template parameter to the SolverCG
// class is the type of the vectors, but the empty angle brackets indicate
// that we simply take the default argument (which is
// <code>Vector@<double@></code>):
SolverCG<> solver(solver_control);
// Now solve the system of equations. The CG solver takes a preconditioner
// as its fourth argument. We don't feel ready to delve into this yet, so we
// tell it to use the identity operation as preconditioner:
solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
// Now that the solver has done its job, the solution variable contains the
// nodal values of the solution function.
}
// @sect4{Step3::output_results}
// The last part of a typical finite element program is to output the results
// and maybe do some postprocessing (for example compute the maximal stress
// values at the boundary, or the average flux across the outflow, etc). We
// have no such postprocessing here, but we would like to write the solution
// to a file.
void Step3::output_results() const
{
// To write the output to a file, we need an object which knows about output
// formats and the like. This is the DataOut class, and we need an object of
// that type:
DataOut<3> data_out;
// Now we have to tell it where to take the values from which it shall
// write. We tell it which DoFHandler object to use, and the solution vector
// (and the name by which the solution variable shall appear in the output
// file). If we had more than one vector which we would like to look at in
// the output (for example right hand sides, errors per cell, etc) we would
// add them as well:
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "solution");
// After the DataOut object knows which data it is to work on, we have to
// tell it to process them into something the back ends can handle. The
// reason is that we have separated the frontend (which knows about how to
// treat DoFHandler objects and data vectors) from the back end (which knows
// many different output formats) and use an intermediate data format to
// transfer data from the front- to the backend. The data is transformed
// into this intermediate format by the following function:
data_out.build_patches();
// Now we have everything in place for the actual output. Just open a file
// and write the data into it, using GNUPLOT format (there are other
// functions which write their data in postscript, AVS, GMV, or some other
// format):
std::ofstream output("solution.vtk");
data_out.write_vtk(output);
}
// @sect4{Step3::run}
// Finally, the last function of this class is the main function which calls
// all the other functions of the <code>Step3</code> class. The order in which
// this is done resembles the order in which most finite element programs
// work. Since the names are mostly self-explanatory, there is not much to
// comment about:
void Step3::run()
{
make_grid();
setup_system();
assemble_system();
solve();
output_results();
}
// @sect3{The <code>main</code> function}
// This is the main function of the program. Since the concept of a
// main function is mostly a remnant from the pre-object oriented era
// before C++ programming, it often does not do much more than
// creating an object of the top-level class and calling its principle
// function.
//
// Finally, the first line of the function is used to enable output of
// some diagnostics that deal.II can generate. The @p deallog
// variable (which stands for deal-log, not de-allog) represents a
// stream to which some parts of the library write output. For
// example, iterative solvers will generate diagnostics (starting
// residual, number of solver steps, final residual) as can be seen
// when running this tutorial program.
//
// The output of @p deallog can be written to the console, to a file,
// or both. Both are disabled by default since over the years we have
// learned that a program should only generate output when a user
// explicitly asks for it. But this can be changed, and to explain how
// this can be done, we need to explain how @p deallog works: When
// individual parts of the library want to log output, they open a
// "context" or "section" into which this output will be placed. At
// the end of the part that wants to write output, one exits this
// section again. Since a function may call another one from within
// the scope where this output section is open, output may in fact be
// nested hierarchically into these sections. The LogStream class of
// which @p deallog is a variable calls each of these sections a
// "prefix" because all output is printed with this prefix at the left
// end of the line, with prefixes separated by colons. There is always
// a default prefix called "DEAL" (a hint at deal.II's history as the
// successor of a previous library called "DEAL" and from which the
// LogStream class is one of the few pieces of code that were taken
// into deal.II).
//
// By default, @p logstream only outputs lines with zero prefixes --
// i.e., all output is disabled because the default "DEAL" prefix is
// always there. But one can set a different maximal number of
// prefixes for lines that should be output to something larger, and
// indeed here we set it to two by calling
// LogStream::depth_console(). This means that for all screen output,
// a context that has pushed one additional prefix beyond the default
// "DEAL" is allowed to print its output to the screen ("console"),
// whereas all further nested sections that would have three or more
// prefixes active would write to @p deallog, but @p deallog does not
// forward this output to the screen. Thus, running this example (or
// looking at the "Results" section), you will see the solver
// statistics prefixed with "DEAL:CG", which is two prefixes. This is
// sufficient for the context of the current program, but you will see
// examples later on (e.g., in step-22) where solvers are nested more
// deeply and where you may get useful information by setting the
// depth even higher.
int main()
{
deallog.depth_console(2);
Step3 laplace_problem;
laplace_problem.run();
return 0;
}'수치해석 > Deal.II' 카테고리의 다른 글
| step3 코드정리 (0) | 2020.11.25 |
|---|---|
| deal.II 주요 class 정리 (0) | 2020.11.18 |
| [VisIt 3.1.3] 설치하기 [visualization tool] [deal ii] [lecture 11] (0) | 2020.11.12 |
| [Deal.ii] [lecture 4] [FEM 구현] (0) | 2020.11.11 |
| lecture3: Deal ll 설치 방법 (3) | 2020.08.17 |