#include <stddef.h>
#include <stdlib.h>
#include <mpi.h>
#include <math.h>
#include "libpme.h"
#include "petscsles.h"
/* fonksiyon prototipleri */
void elem_tri(double *x,double *y,double *ek[3],double ef[3]) ;
int fem_poisson(pMesh mesh,int mypid) ;
double f(double x,double y) ;
double exact(double x,double y) ;
int refine_uniform(pMesh) ;
main(
int argc,
char *argv[])
{
/* tanimlar */
pMesh mesh ;
char filename[256] ;
int rc,i ;
int mypid, numprocs ;
double startwtime ;
double endwtime ;
/* MPI ile PETSC'yi ilk kullanima hazirla */
MPI_Init(&argc,&argv) ;
MPI_Comm_rank(MPI_COMM_WORLD,&mypid) ;
MPI_Comm_size(MPI_COMM_WORLD,&numprocs) ;
PetscInitialize(&argc,&argv,(char *)0,0) ;
/* 1.**** Grid'i oku *****************************************/
sprintf(filename,"%s.%d.msh",argv[1],mypid) ;
mesh = pme_new_msh(MPI_COMM_WORLD,NULL) ;
rc = pme_read_msh(mesh,filename) ;
/* Zamanlamayi baslat */
MPI_Barrier(MPI_COMM_WORLD) ;
startwtime = MPI_Wtime();
/* 2.**** Grid'i parcala ************************************/
pme_part_inertia(mesh) ;
/* 3. *** Grid'e ince uyarlama yap ***************************/
for(i=0 ; i < atoi(argv[2]) ; i++) refine_uniform(mesh) ;
printf("%d\n",pme_num_entity(mesh,PME_FACE)) ;
/* 4 *** Sonlu Eleman yontemi ile Laplace denklemini coz ****/
fem_poisson(mesh,mypid) ;
/* Zamanlamayi bitir */
endwtime = MPI_Wtime();
printf("Zaman = %f\n",endwtime-startwtime) ;
MPI_Finalize() ;
exit(0) ;
}
int fem_poisson(
pMesh mesh,
int mypid)
{
void *temp1 ;
void *temp2 ;
pFace f ;
pVertex v ;
int ids[3] ;
int i,j ;
double x[3] ;
double y[3] ;
double ef[3] ; /* eleman sag taraf */
double arry[9] ; /* eleman matris */
double *ek[3] ; /* eleman matris 2d arry referansi */
pEnList enlist ;
Vec U, b, u ; /* hesaplanmis cozum, sag taraf,
gercek cozum */
Mat A ; /* dogrusal sistem matrisi */
SLES sles ; /* SLES ortami */
int ierr ;
int n,N ; /* sahip olunan yorel varlik sayisi
global varlik sayisi */
int its ; /* iterasyon sayisi */
double norm ;
Scalar neg_one = -1.0,
zero=0.0,
one=1.0 ;
IS is ;
int *rows ;
int n_gb ;
double *gb_f ;
Scalar *solnu ;
Scalar *solnU ;
double *d ;
int *indxs ;
pEntity owner_entity ;
int owner_pid ;
pFace face ;
pVertex vertex ;
ek[0] = &(arry[0]) ; ek[1] = &(arry[3]) ; ek[2] = &(arry[6]) ;
/* 1 ** Vektorleri yarat ******************************/
pme_glob_numbering(mesh,PME_VERTEX) ;
pme_glob_num_entity(mesh,PME_VERTEX,&n,&N) ;
MYMALLOC(indxs,int *,n*sizeof(int)) ;
MYMALLOC(solnu,Scalar *,n*sizeof(Scalar)) ;
ierr = VecCreateMPI(MPI_COMM_WORLD,n,N,&U) ;
ierr = VecDuplicate(U,&b) ;
ierr = VecDuplicate(U,&u) ;
ierr = VecSet(&zero,U) ;
/* 2 *** Matrisi yarat, eleman hesaplari yap,
ve matris kurulumunu yap *******/
ierr = MatCreateMPIAIJ(MPI_COMM_WORLD,n,n,N,N,
0,PETSC_NULL,0,PETSC_NULL,&A);
/* Elemen hesaplari */
temp1 = NULL ;
while(face = pme_next_entity(mesh,PME_FACE,&temp1)) {
i = 2 ;
temp2 = NULL ;
pme_get_reln(mesh,face,PME_VERTEX,&enlist) ;
while(vertex = pme_next_enlist(enlist,&temp2) ) {
ids[i] = pme_glob_number(mesh,vertex) - 1 ;
x[i] = pme_get_vertex_xyz(vertex,0) ;
y[i] = pme_get_vertex_xyz(vertex,1) ;
i-- ;
}
pme_free_enlist(enlist) ;
elem_tri(x,y,ek,ef) ;
ierr = MatSetValues(A,3,ids,3,ids,arry,ADD_VALUES);
VecSetValues(b,3,ids,ef,ADD_VALUES);
}
/* matris kurulumu */
ierr = VecAssemblyBegin(b) ; CHKERRQ(ierr);
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = VecAssemblyEnd(b) ; CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY) ; CHKERRQ(ierr);
/** 3 * Dirichlet sinir kosullari *********/
n_gb = pme_num_gb_entity(mesh,PME_VERTEX) ;
MYMALLOC(rows,int *,n_gb*sizeof(int)) ;
MYMALLOC(gb_f,double *,n_gb*sizeof(double)) ;
temp1 = NULL ;
i = 0 ;
while(v = pme_next_gb_entity(mesh,PME_VERTEX,&temp1) ) {
pme_get_pb_owner(mesh,v,&owner_pid,&owner_entity) ;
if (owner_pid == mypid) {
rows[i] = pme_glob_number(mesh,v) - 1 ;
gb_f[i] = exact(pme_get_vertex_xyz(v,0),
pme_get_vertex_xyz(v,1) ) ;
i++ ;
}
}
n_gb = i ;
VecSetValues(b,n_gb,rows,gb_f,INSERT_VALUES);
ierr = ISCreateGeneral(MPI_COMM_SELF,n_gb,rows,&is); CHKERRQ(ierr);
ierr = MatZeroRows(A,is,&one); CHKERRQ(ierr);
ierr = ISDestroy(is); CHKERRQ(ierr);
/** 4 * Dogrusal denklemi coz ************************/
ierr = SLESCreate(MPI_COMM_WORLD,&sles) ; CHKERRQ(ierr);
ierr = SLESSetOperators(sles,A,A,
DIFFERENT_NONZERO_PATTERN); CHKERRQ(ierr);
ierr = SLESSetFromOptions(sles) ; CHKERRQ(ierr);
ierr = SLESSolve(sles,b,U,&its) ; CHKERRQ(ierr);
/** 5 * Hata'yi kontrol et ***********************/
temp1 = NULL ;
i = 0 ;
while(v = pme_next_entity(mesh,PME_VERTEX,&temp1)) {
pme_get_pb_owner(mesh,v,&owner_pid,&owner_entity);
if(owner_pid == mypid) {
solnu[i] = exact(pme_get_vertex_xyz(v,0),
pme_get_vertex_xyz(v,1) );
indxs[i] = pme_glob_number(mesh,v) - 1;
i++ ;
}
}
n = i ;
VecSetValues(u,n,indxs,solnu,INSERT_VALUES);
ierr = VecAssemblyBegin(u); CHKERRQ(ierr);
ierr = VecAssemblyEnd(u) ; CHKERRQ(ierr);
ierr = VecAXPY(&neg_one,u,U); CHKERRQ(ierr);
ierr = VecNorm(U,NORM_MAX,&norm);
printf("Norm= %lE Iter=%d\n",norm,its) ;
ierr = VecDestroy(U);
ierr = VecDestroy(u);
ierr = VecDestroy(b);
ierr = MatDestroy(A);
ierr = SLESDestroy(sles);
}
/*
* Dogrusal sekil (shape) fonksiyonlari ile ucgen uzerinde
* eleman hesaplarini yap
*/
void elem_tri(
double *x,
double *y,
double *ek[3],
double ef[3])
{
static double psi[3] = {0.33333333333,0.33333333333,0.33333333333} ;
static double weight = 0.5 ;
double x_l, y_l ;
double detj, one_detj,fac ;
int i,j,k ;
double p1, p2 ;
int k1, k2 ;
/* x,y */
x_l = psi[0]*x[0] + psi[1]*x[1] + psi[2]*x[2] ;
y_l = psi[0]*y[0] + psi[1]*y[1] + psi[2]*y[2] ;
/* Jakobian */
detj = (x[1]-x[0])*(y[2]-y[0]) - (x[2]-x[0])*(y[1]-y[0]) ;
one_detj = 1.0/detj ;
/* Integral hesaplamasi */
for(i=0 ; i < 3 ; i++) {
ef[i] = f(x_l,y_l)*psi[i]*detj*weight ;
k1 = (i+1)%3 ; k2 = (i+2)%3 ;
p1 = (y[k1]-y[k2]) ;
p2 = (x[k1]-x[k2]) ;
ek[i][i] = -weight*one_detj*( p1*p1 + p2*p2 ) ;
for(j=i+1 ; j < 3 ; j++) {
k = 3 - i - j ;
ek[j][i] =
ek[i][j] = -weight*one_detj*
((y[j]-y[k])*(y[k]-y[i])+(x[k]-x[j])*(x[i]-x[k])) ;
}
}
}
double f(
double x,
double y)
{
double t ;
return(0.0) ;
}
double exact(
double x,
double y)
{
return(x*y) ;
}
/* Bu fonksiyon butun elemenlara ince uyarlama yapar */
refine_uniform(
pMesh mesh)
{
AttchTag mark_tag ;
pFace face ;
void *temp ;
pme_new_attch_tag(mesh,NULL,PME_ATTCH_NOT_FREE,&mark_tag) ;
/* butun elemanlari ince uyarlama icin isaretle */
temp = NULL ;
while (face = pme_next_entity(mesh,PME_FACE,&temp) ) {
pme_ins_data_attch(face,mark_tag,(void *)1) ;
}
/* Paralel Rivara ince uyarlama fonksiyonunu cagir */
pme_ref_rivara(mesh,mark_tag) ;
/* isaretleri sil */
temp = NULL ;
while( face = pme_next_entity(mesh,PME_FACE,&temp) ) {
pme_del_data_attch(face,mark_tag) ;
}
pme_del_attch_tag(mesh,mark_tag) ;
}