1. 有限元方法大家可以找一本书普及一个基本知识。
2. 写一个有限元方法的程序可以参考
/home/tlu/Downloads/trilinos-10.6.1/packages/didasko/examples/epetra/ex13.cpp
这个程序没有处理边界条件， 我们这儿的边界条件是Dirichlet 边界， 边界上的函数值等于0.
因此，只要把生成矩阵的对应边界的行对角元置成1， 其他元素置成零即可。
3. 得到的线性方程组的求解可以参考
http://hi.baidu.com/motioo/blog/item/4b933a2719795601908f9d49.htm
4. Ml preconditioner 可以参考
/home/tlu/Downloads/trilinos-10.6.1/packages/didasko/examples/ml/ex1.cpp
5. 要写这个程序还要easymesh 生成网格, 请参考
http://www-dinma.univ.trieste.it/nirftc/research/easymesh/
需要使用easymesh 生成的三角形的节点的坐标和三角形的三个节点的编号。
6 计算结果的画图，
3 -------------4
2
0 --------------1
5点， 4个三角形， 分别是
0 1 2 ;
1 2 4 ;
4 2 3 ;
3 2 0.
针对上面的剖分，假设四个角点的函数值是0，中间点的函数值是10
则对应tecplot的文件格式举例如下
a.tec
文件开始
TITLE = "An example show to how plot data by using tecplot"
VARIABLES = "X","Y", "u"
ZONE T="SQUARE",N=5, F=FEPOINT, ET=QUADRILATERAL
0 0 0
0 1 0
0.5 0.5 10
1 1 0
0 1 0
1 2 3
2 3 5
5 3 4
4 3 0
a.tec 结束，
注意： 三角形中的顶点的编号是从1 开始的。
7. 原来写的两个程序， 一个是把 easymesh 输出的网格转成 .msh
，另一个是把 .msh 转成 tecplot 文件。
#include <iostream>
#include <fstream>
#include <string>
#include <assert.h>
/* you have easymesh a.d, a.n a.s a.e,
* ./easy2msh a a.msh
* and you can 2d msh file a.msh
*/
int main(int argc, char ** argv){
int n, i,j ;
std::string name , filename;
std::string outfile;
std::ifstream is;
char ch;
double x;
if(argc < 2) {
std::cout<<"input in fileanme "<<std::endl;
std::cin >> name;
}
else{
name = argv[1];
}
if(argc < 3) {
std::cout<<"input out fileanme "<<std::endl;
std::cin >> outfile;
}
else{
outfile = argv[2];
}
filename = name + ".n";
is.open(filename.c_str());
assert(is);
is >> n;
std::ofstream os(outfile.c_str());
os.precision(15);
os.setf(std::ios::scientific);
os<<"Dimension"<<std::endl;
os<<2<<std::endl;
os<<"Vertices"<<std::endl;
os<<n<<std::endl;
for( i = 0; i < n; i ++) {
is >> j;
is >> ch;
is >> x;
os << x <<" ";
is >> x;
os << x <<" ";
is >> j;
os << j <<std::endl;
}
/*
is.close();
filename = name +".s";
is.open(filename.c_str());
is >> n;
os<<"Edges"<<std::endl;
os<<n<<std::endl;
for( i = 0; i < n; i ++){
is >> j;
is >> ch;
is >>j; j ++;
os <<j << " ";
is >>j; j ++;
os <<j <<" ";
is >>j;
is >>j;
is >>j; j ++;
is >>j;
os <<j<<std::endl;
}
*/
is.close();
filename = name +".e";
is.open(filename.c_str());
is >> n;
os <<"Triangles"<<std::endl;
os<<n<<std::endl;
for( i = 0; i < n; i ++) {
is >> j;
is >> ch;
is >> j; j++;
os << j <<" ";
is >> j;j++;
os << j <<" ";
is >> j; j++;
os << j <<" ";
is >> j;is >> j;is >> j;is >> j;is >> j;is >> j;
is >> x ; is >> x;
is >> j ;
os << j <<std::endl;
}
is.close();
os <<"End"<<std::endl;
os.close();
return 1;
}
/* Author: Lu Tiao Email: tlu@math.pku.edu.cn */
/* Translate the msh file to tec file */
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <string>
int main(int argc, char * argv[]){
if( argc <3) {
std::cout<<"1.input the msh file "<<std::endl;
std::cout<<"2. intput the tec file " <<std::endl;
exit(1);
}
int i,j,k, n_points, n_triangles, n_quadriaterals;
std::string dummy;
double coord;
std::ifstream is(argv[1]);
std::ofstream os(argv[2]);
is.precision(15);
os.precision(15);
os<<"TITLE = \"Mesh\""<<std::endl;
os<<"VARIABLES = \"X\",\"Y\""<<std::endl;
os<<"ZONE T=\"SQUARE\",N=";
do{
getline(is,dummy);
if(dummy == "Vertices" ){
is >> n_points;
os<<n_points<<", F=FEPOINT, ET=QUADRILATERAL"<<std::endl;
for( i = 0; i < n_points; i ++) {
for( j = 0; j < 2; j ++) {
is >> coord;
os<< coord<<" ";
}
os <<std::endl;
is >> k;
}
}
if( dummy == "Triangles" ){
is >> n_triangles;
for( i = 0; i < n_triangles; i ++) {
for( k = 0; k < 3; k ++) {
is >> j;
os << j <<" ";
}
os << j << std::endl;
is >>k ;
}
}
if( dummy == "Quadrilaterals" ){
is >> n_quadriaterals;
for( i = 0; i < n_quadriaterals; i ++) {
for( k = 0; k < 4; k ++) {
is >> j;
os << j <<" ";
}
os << std::endl;
is >>k ;
}
}
}while(!is.eof());
is.close();
os.close();
return 1;
}
