Index: configure.ac =================================================================== --- configure.ac (revision 5e1467a198e04b81d706d65ad5db6f21710e8e06) +++ configure.ac (revision 4bce0762b411f260032c6e4d1f6bf6f607e6ef29) @@ -166,5 +166,6 @@ src/Makefile src/main/Makefile - src/ext/Makefile + src/ext/C/Makefile + src/ext/triangle/Makefile doc/Doxyfile doc/Makefile Index: src/datatypes/Mesh.F90 =================================================================== --- src/datatypes/Mesh.F90 (revision 6ca9419e7f944c485368f2f83012dc0854db08b6) +++ src/datatypes/Mesh.F90 (revision 438cb4851fe5a0091c37a009b8d5f75554f6182d) @@ -243,52 +243,4 @@ end function Mesh_newInit - - !---------------------------------------------------------- - !! Allocates Mesh's main data - !---------------------------------------------------------- - subroutine Mesh_allocate(M, & - nfrelt, & - mhead, & - freemap, & - coords, & - eptnmap, & - freemask) - implicit none - - type(Mesh), intent(in out) :: M - logical, intent(in), optional :: nfrelt, mhead, freemap - logical, intent(in), optional :: coords, eptnmap, freemask - logical :: all=.false. - if (.not.present(nfrelt) .and.& - (.not.present(mhead)) .and.& - (.not.present(freemap)).and.& - (.not.present(coords)) .and.& - (.not.present(eptnmap)) .and.& - (.not.present(freemask))) then - if (.not.associated(M%nfrelt)) allocate(M%nfrelt(M%nell)) - if (.not.associated(M%mhead)) allocate(M%mhead(M%mfrelt,M%nell)) - if (.not.associated(M%freemap)) allocate(M%freemap(M%ngf)) - if (.not.associated(M%coords)) allocate(M%coords(M%nsd,M%nnode)) - if (.not.associated(M%eptnmap)) allocate(M%eptnmap(M%nell)) - if (.not.associated(M%freemask)) allocate(M%freemask(M%ngf)) - all = .true. - end if - - if (present(nfrelt) .and.(.not.all).and.& - (.not.associated(M%nfrelt))) allocate(M%nfrelt(M%nell)) - if (present(mhead) .and.(.not.all).and.& - (.not.associated(M%mhead))) allocate(M%mhead(M%mfrelt,M%nell)) - if (present(freemap) .and.(.not.all).and.& - (.not.associated(M%freemap))) allocate(M%freemap(M%ngf)) - if (present(coords) .and.(.not.all).and.& - (.not.associated(M%coords))) allocate(M%coords(M%nsd,M%nnode)) - if (present(eptnmap) .and.(.not.all).and.& - (.not.associated(M%eptnmap))) allocate(M%eptnmap(M%nell)) - if (present(freemask).and.(.not.all).and.& - (.not.associated(M%freemask))) allocate(M%freemask(M%ngf)) - - end subroutine Mesh_allocate - - !------------------------- !! Destructor @@ -387,5 +339,5 @@ nnode=ngf call Mesh_Init(M, nell, ngf, nsd, mfrelt, nnode) - call Mesh_allocate(M,coords=.true.,freemap=.true.) ! needing the coords... + allocate(M%coords(M%nsd,M%nnode),M%freemap(M%ngf)) ! needing the coords... M%freemap= (/ (i,i=1,ngf) /) do j=1,n @@ -1503,4 +1455,5 @@ end do else + if (.NOT.associated(M%nfrelt)) allocate(M%nfrelt(M%nell)) call MPI_RECV(M%nfrelt, M%nell, MPI_INTEGER, & D_MASTER, D_TAG_MESH_NFRELT, MPI_COMM_WORLD, status, ierr) @@ -1516,4 +1469,5 @@ end do else + if (.NOT.associated(M%mhead)) allocate(M%mhead(M%mfrelt,M%nell)) call MPI_RECV(M%mhead, M%mfrelt*M%nell, MPI_INTEGER, & D_MASTER, D_TAG_MESH_MHEAD, MPI_COMM_WORLD, status, ierr) @@ -1529,4 +1483,5 @@ end do else + if (.not.associated(M%freemap)) allocate(M%freemap(M%ngf)) call MPI_RECV(M%freemap, M%ngf, MPI_INTEGER, & D_MASTER, D_TAG_MESH_FREEMASK, MPI_COMM_WORLD, status, ierr) @@ -1542,4 +1497,5 @@ end do else + if (.not.associated(M%coords)) allocate(M%coords(M%nsd,M%nnode)) call MPI_RECV(M%coords, M%nnode*M%nsd, MPI_xyzkind, & D_MASTER, D_TAG_MESH_FREEMASK, MPI_COMM_WORLD, status, ierr) @@ -1549,4 +1505,7 @@ ! Send/recv 'eptnmap' if (present(eptnmap)) then + call MPI_BCAST(M%nparts, 1, MPI_INTEGER, & + D_MASTER, MPI_COMM_WORLD, ierr) + if (ismaster()) then do p = 1,numprocs-1 @@ -1555,6 +1514,8 @@ end do else + if (.not.associated(M%eptnmap)) allocate(M%eptnmap(M%nell)) call MPI_RECV(M%eptnmap, M%nell, MPI_INTEGER, & D_MASTER, D_TAG_MESH_EPTNMAP, MPI_COMM_WORLD, status, ierr) + M%parted = .true. end if end if @@ -1592,4 +1553,5 @@ end do else + if (.not.associated(M%freemask)) allocate(M%freemask(M%ngf)) call MPI_RECV(M%freemask, M%ngf, MPI_BYTE, & D_MASTER, D_TAG_MESH_FREEMASK, MPI_COMM_WORLD, status, ierr) @@ -1685,4 +1647,33 @@ !! Plotting routines ! + + !> Show a sequence of plots of the mesh. + subroutine Mesh_pl2D_mesh(Msh) + type(Mesh), intent(inout) :: Msh + + call Mesh_pl2D_pointCloud(Msh,D_PLPLOT_INIT) + ! Plots mesh's dual graph + call Mesh_pl2D_plotGraphDual(Msh,D_PLPLOT_END) + ! Mesh & its Dual Graph + call Mesh_pl2D_plotMesh(Msh, D_PLPLOT_INIT) + call Mesh_pl2D_plotGraphDual(Msh, D_PLPLOT_END) + end subroutine Mesh_pl2D_mesh + + !> Show plots of the mesh partitions. + subroutine Mesh_pl2D_partitions(Msh) + type(Mesh), intent(inout) :: Msh + + ! Draw colored partitoined graph + call Mesh_pl2D_plotGraphParted(Msh) + + ! Plot partitions of the mesh + ! NB: Check for multivariable case! TODO + call Mesh_pl2D_Partition(Msh) + ! Partition with Dual Graph upon it + call Mesh_pl2D_Partition(Msh, D_PLPLOT_INIT) + call Mesh_pl2D_plotGraphDual(Msh, D_PLPLOT_CONT) + call Mesh_pl2D_pointCloud(Msh,D_PLPLOT_END) + end subroutine Mesh_pl2D_partitions + !------------------------------------------------ !! Plot cloud field @@ -2524,5 +2515,8 @@ end do - if (M%nfrelt(e) == 1) then ! 1-node boundary element + if (M%nfrelt(e) == 0) then + cycle + + else if (M%nfrelt(e) == 1) then ! 1-node boundary element ! Fake centre of 1-node element by simply assigning to Index: src/ext/C/Makefile.am =================================================================== --- src/ext/C/Makefile.am (revision dd63f9fdf01c73878366fca9ec399ccab07b12cc) +++ src/ext/C/Makefile.am (revision dd63f9fdf01c73878366fca9ec399ccab07b12cc) @@ -0,0 +1,15 @@ +lib_LTLIBRARIES = libdoug-ext.la +libdoug_ext_la_SOURCES = doug_ext.f90 ops.c +libdoug_ext_la_FCFLAGS=-I@top_builddir@/src -I@top_builddir@/src/main +noinst_PROGRAMS=c_doug_aggr c_doug_cg c_doug_pcg1 + +c_doug_aggr_SOURCES=main_aggr.c +# this file is compiled with 'mpicc' and in static linking does not find MPI library for Fortran, +# thus we hack it by providing required libraries, however it only works for OpenMPI +c_doug_aggr_LDADD=libdoug-ext.la @top_builddir@/src/libdoug.la @top_builddir@/src/main/doug_aggr-main_drivers.o -lmpi_f90 -lmpi_f77 + +c_doug_cg_SOURCES=main_cg.c +c_doug_cg_LDADD=libdoug-ext.la @top_builddir@/src/libdoug.la @top_builddir@/src/main/doug_aggr-main_drivers.o -lmpi_f90 -lmpi_f77 + +c_doug_pcg1_SOURCES=main_pcg1.c +c_doug_pcg1_LDADD=libdoug-ext.la @top_builddir@/src/libdoug.la @top_builddir@/src/main/doug_aggr-main_drivers.o -lmpi_f90 -lmpi_f77 Index: src/ext/C/doug_ext.f90 =================================================================== --- src/ext/C/doug_ext.f90 (revision dd63f9fdf01c73878366fca9ec399ccab07b12cc) +++ src/ext/C/doug_ext.f90 (revision dd63f9fdf01c73878366fca9ec399ccab07b12cc) @@ -0,0 +1,147 @@ +!#include + +!> Initialize DOUG, init_type: 1 - parallel, 2 - serial. +subroutine ext_DOUG_Init(init_type) + use DOUG_utils + integer, intent(in), optional :: init_type + call DOUG_Init(init_type) +end subroutine ext_DOUG_Init + +subroutine ext_DOUG_Finalize() + use DOUG_utils + call DOUG_Finalize() +end subroutine ext_DOUG_Finalize + +! Allocate memory for SpMtx class. +subroutine alloc_SpMtx(A) + use SpMtx_class + type(SpMtx),pointer :: A + allocate(A) + A = SpMtx_New() +end subroutine alloc_SpMtx + +! Allocate memory for Mesh class. +subroutine alloc_Mesh(M) + use Mesh_class + type(Mesh),pointer :: M + allocate(M) + M = Mesh_New() +end subroutine alloc_Mesh + +subroutine ext_parallelDistributeInput(M, A, b, nparts, A_ghost, nlf) + use main_drivers + use globals, only: sctls + implicit none + type(Mesh), intent(in out) :: M + type(SpMtx), intent(out) :: A + real(kind=rk), intent(out) :: b(*) + integer, intent(in) :: nparts + type(SpMtx),intent(in out) :: A_ghost + integer, intent(out) :: nlf + + real(kind=rk), dimension(:), pointer :: b_ + integer, dimension(6) :: part_opts + + call parallelDistributeInput(sctls%input_type, M, A, b_, nparts, part_opts, A_ghost) + b(:size(b_)) = b_ + nlf = size(b_) + deallocate(b_) + +end subroutine ext_parallelDistributeInput + +subroutine ext_cg(A, b, M, xl, nlf) + use cg_mod + implicit none + + type(SpMtx),intent(in out) :: A ! System matrix (sparse) + real(kind=rk),dimension(nlf) :: b ! RHS + real(kind=rk),dimension(nlf) :: xl ! Solution + ! Mesh - aux data for Ax operation + type(Mesh),intent(in) :: M + integer, intent(in) :: nlf + + ! Solution statistics + type(ConvInf) :: solinf + + real(kind=rk),pointer :: b_(:), xl_(:) ! RHS + allocate(b_(size(b))) + b_ = b + allocate(xl_(size(xl))) + xl_ = xl + + call cg(A, b_, xl_, M, solinf=solinf) + + xl = xl_ + deallocate(b_) + deallocate(xl_) +end subroutine ext_cg + +subroutine ext_pmvmCommStructs_init(A,M) + use SpMtx_operation + + type(SpMtx), intent(in) :: A + type(Mesh), intent(in) :: M + + call pmvmCommStructs_init(A,M) +end subroutine ext_pmvmCommStructs_init + +subroutine ext_SpMtx_pmvm(r,A,x,M,n) + use SpMtx_operation + + type(SpMtx), intent(in) :: A + type(Mesh), intent(in) :: M + integer, intent(in) :: n + real(kind=rk), intent(in), target :: x(n) + real(kind=rk), intent(out), target :: r(n) + + real(kind=rk), pointer :: x_(:), r_(:) + + x_ => x + r_ => r + call SpMtx_pmvm(r_,A,x_,M) + +end subroutine ext_SpMtx_pmvm + +subroutine ext_vect_dot(v,x,y,n) + use Vect_mod + + real(kind=rk), intent(in), target :: x(n), y(n) + integer, intent(in) :: n + real(kind=rk), intent(out) :: v + + real(kind=rk), pointer :: x_(:), y_(:) + + x_ => x + y_ => y + v = Vect_dot_product(x_,y_) + +end subroutine ext_vect_dot + +subroutine ext_preconditioner_1level(sol,A,rhs,M,A_interf_,refactor_,nlf) + use pcg_mod + + implicit none + real(kind=rk),dimension(nlf),target :: sol !< solution + type(SpMtx) :: A !< sparse system matrix + real(kind=rk),dimension(nlf),target :: rhs !< right hand side + type(Mesh),intent(in) :: M !< Mesh + type(SpMtx),optional :: A_interf_ !< matr@interf. + logical,intent(inout) :: refactor_ + integer, intent(in) :: nlf + + !< residual vector, allocated + !! here for multiplicative Schwarz + real(kind=rk),dimension(:),pointer :: res => NULL() + + real(kind=rk), pointer :: sol_(:), rhs_(:) + + sol_ => sol + rhs_ => rhs + sol_ = 0 + call preconditioner_1level(sol_,A,rhs_,M,res,A_interf_,refactor_) + + if (sctls%method/=0) then + call Add_common_interf(sol_,A,M) + endif + +end subroutine ext_preconditioner_1level Index: src/ext/C/doug_ext.h =================================================================== --- src/ext/C/doug_ext.h (revision dd63f9fdf01c73878366fca9ec399ccab07b12cc) +++ src/ext/C/doug_ext.h (revision dd63f9fdf01c73878366fca9ec399ccab07b12cc) @@ -0,0 +1,15 @@ +#ifndef DOUG_EXT_H +#define DOUG_EXT_H + +extern void ext_doug_init_(int *init_type); +extern void ext_doug_finalize_(); +extern void alloc_spmtx_(void **A); +extern void alloc_mesh_(void **M); +extern void ext_paralleldistributeinput_(void *M, void *A, double *b, int *nparts, void *A_ghost, int *nlf); +extern void ext_cg_(void *A, double *b, void *M, double *xl, int *nlf); +extern void ext_spmtx_pmvm_(double *b, void *A, double *xl, void *M, int *nlf); +extern double ext_vec_dot_(double *v, double *x, double *y, int *nlf); +extern void ext_pmvmcommstructs_init_(void *A, void *M); +extern void ext_preconditioner_1level_(double *sol, void *A, double *rhs, void *M, void *A_interf_, int *refactor_, int *nlf); + +#endif Index: src/ext/C/main_aggr.c =================================================================== --- src/ext/C/main_aggr.c (revision dd63f9fdf01c73878366fca9ec399ccab07b12cc) +++ src/ext/C/main_aggr.c (revision dd63f9fdf01c73878366fca9ec399ccab07b12cc) @@ -0,0 +1,37 @@ +#include +#include + +#include "doug_ext.h" + +int main(int argc, char **argv) +{ + int i; + void *A, *A_ghost; + void *M; + double *xl; + double *b; + int nparts=1; + int init_type=1; + int nlf; + + xl = malloc(225*sizeof(double)); + b = malloc(225*sizeof(double)); + + ext_doug_init_(&init_type); // parallel + alloc_spmtx_(&A); + alloc_spmtx_(&A_ghost); + alloc_mesh_(&M); + printf("Initialized, matrix %x, mesh %x\n", A, M); + ext_paralleldistributeinput_(M, A, b, &nparts, A_ghost, &nlf); + + ext_cg_(A, b, M, xl, &nlf); + + ext_doug_finalize_(); + + printf("result size=%d:\n", nlf); + for(i=0; i +#include +#include + +#include + +#include "doug_ext.h" +#include "ops.h" + +static void cg(void *A, void *M, double *xl, double *b, int nlf) +{ + double *r = malloc(nlf*sizeof(double)); + double *p = malloc(nlf*sizeof(double)); + double *q = malloc(nlf*sizeof(double)); + double tol, rho, rho_old, alpha, beta, tmp; + int it; + int rank, i; + + tol = 1E-12; + + // initialize pmvm constructs for communication + ext_pmvmcommstructs_init_(A, M); + + rho_old = 1.0; + ext_spmtx_pmvm_(r,A,xl,M, &nlf); + daxpy(r, -1.0, r, b, nlf); + for(i=0; i tol*tol) { + ext_vect_dot_(&rho, r, r, &nlf); + if (rank==0) + printf("rank=%d, it = %d, sqrt(rho) = %e\n", rank, it, sqrt(rho)); + + beta = rho/rho_old; + daxpy(p, beta, p, r, nlf); + ext_spmtx_pmvm_(q,A,p,M, &nlf); + ext_vect_dot_(&tmp, p,q,&nlf); + alpha = rho/tmp; + daxpy(xl, alpha, p, xl, nlf); + daxpy(r, -alpha, q, r, nlf); + rho_old = rho; + it++; + } +} + +int main(int argc, char **argv) +{ + int i; + void *A, *A_ghost; + void *M; + double *xl; + double *b; + int nparts=1; + int init_type=1; + int nlf; + + xl = malloc(225*sizeof(double)); + b = malloc(225*sizeof(double)); + + ext_doug_init_(&init_type); // parallel + alloc_spmtx_(&A); + alloc_spmtx_(&A_ghost); + alloc_mesh_(&M); + printf("Initialized, matrix %x, mesh %x\n", A, M); + ext_paralleldistributeinput_(M, A, b, &nparts, A_ghost, &nlf); + + cg(A, M, xl, b, nlf); + + ext_doug_finalize_(); + + for(i=0; i +#include +#include + +#include + +#include "doug_ext.h" +#include "ops.h" + +static void pcg(void *A, void *M, void *A_ghost, double *xl, double *b, int nlf) +{ + double *r = malloc(nlf*sizeof(double)); + double *p = malloc(nlf*sizeof(double)); + double *q = malloc(nlf*sizeof(double)); + double *z = malloc(nlf*sizeof(double)); + double tol, rho, rho_old, alpha, beta, tmp; + double ratio_norm, init_norm, res_norm; + int it, refactor=1; + int rank, i; + + tol = 1E-12; + + // initialize pmvm constructs for communication + ext_pmvmcommstructs_init_(A, M); + + rho_old = 1.0; + ext_spmtx_pmvm_(r,A,xl,M, &nlf); + daxpy(r, -1.0, r, b, nlf); + + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + ext_vect_dot_(&init_norm, r,r,&nlf); + if (init_norm == 0.0) init_norm = 1.0; + ratio_norm = 1.0; + + it = 0; + while (ratio_norm > tol*tol) { + it++; + ext_preconditioner_1level_(z,A,r,M,A_ghost,&refactor,&nlf); + refactor = 0; + + ext_vect_dot_(&rho, r, z, &nlf); + + if (it==1) { + ops_copy(p, z, nlf); + } else { + beta = rho/rho_old; + daxpy(p, beta, p, z, nlf); + } + ext_spmtx_pmvm_(q,A,p,M, &nlf); + ext_vect_dot_(&tmp, p,q,&nlf); + alpha = rho/tmp; + daxpy(xl, alpha, p, xl, nlf); + daxpy(r, -alpha, q, r, nlf); + rho_old = rho; + + ext_vect_dot_(&res_norm, r, r, &nlf); + ratio_norm = res_norm / init_norm; + if (rank==0) + printf("rank=%d, it = %d, sqrt(res_norm) = %e\n", rank, it, sqrt(res_norm)); + } +} + +int main(int argc, char **argv) +{ + int i; + void *A, *A_ghost; + void *M; + double *xl; + double *b; + int nparts=1; + int init_type=1; + int nlf; + + xl = malloc(225*sizeof(double)); + b = malloc(225*sizeof(double)); + + ext_doug_init_(&init_type); // parallel + alloc_spmtx_(&A); + alloc_spmtx_(&A_ghost); + alloc_mesh_(&M); + printf("Initialized, matrix %x, mesh %x\n", A, M); + ext_paralleldistributeinput_(M, A, b, &nparts, A_ghost, &nlf); + + pcg(A, M, A_ghost, xl, b, nlf); + + ext_doug_finalize_(); + + for(i=0; i - -!> Initialize DOUG, init_type: 1 - parallel, 2 - serial. -subroutine ext_DOUG_Init(init_type) - use DOUG_utils - integer, intent(in), optional :: init_type - call DOUG_Init(init_type) -end subroutine ext_DOUG_Init - -subroutine ext_DOUG_Finalize() - use DOUG_utils - call DOUG_Finalize() -end subroutine ext_DOUG_Finalize - -! Allocate memory for SpMtx class. -subroutine alloc_SpMtx(A) - use SpMtx_class - type(SpMtx),pointer :: A - allocate(A) - A = SpMtx_New() -end subroutine alloc_SpMtx - -! Allocate memory for Mesh class. -subroutine alloc_Mesh(M) - use Mesh_class - type(Mesh),pointer :: M - allocate(M) - M = Mesh_New() -end subroutine alloc_Mesh - -subroutine ext_parallelDistributeInput(M, A, b, nparts, A_ghost, nlf) - use main_drivers - use globals, only: sctls - implicit none - type(Mesh), intent(in out) :: M - type(SpMtx), intent(out) :: A - real(kind=rk), intent(out) :: b(*) - integer, intent(in) :: nparts - type(SpMtx),intent(in out) :: A_ghost - integer, intent(out) :: nlf - - real(kind=rk), dimension(:), pointer :: b_ - integer, dimension(6) :: part_opts - - call parallelDistributeInput(sctls%input_type, M, A, b_, nparts, part_opts, A_ghost) - b(:size(b_)) = b_ - nlf = size(b_) - deallocate(b_) - -end subroutine ext_parallelDistributeInput - -subroutine ext_cg(A, b, M, xl, nlf) - use cg_mod - implicit none - - type(SpMtx),intent(in out) :: A ! System matrix (sparse) - real(kind=rk),dimension(nlf) :: b ! RHS - real(kind=rk),dimension(nlf) :: xl ! Solution - ! Mesh - aux data for Ax operation - type(Mesh),intent(in) :: M - integer, intent(in) :: nlf - - ! Solution statistics - type(ConvInf) :: solinf - - real(kind=rk),pointer :: b_(:), xl_(:) ! RHS - allocate(b_(size(b))) - b_ = b - allocate(xl_(size(xl))) - xl_ = xl - - call cg(A, b_, xl_, M, solinf=solinf) - - xl = xl_ - deallocate(b_) - deallocate(xl_) -end subroutine ext_cg - -subroutine ext_pmvmCommStructs_init(A,M) - use SpMtx_operation - - type(SpMtx), intent(in) :: A - type(Mesh), intent(in) :: M - - call pmvmCommStructs_init(A,M) -end subroutine ext_pmvmCommStructs_init - -subroutine ext_SpMtx_pmvm(r,A,x,M,n) - use SpMtx_operation - - type(SpMtx), intent(in) :: A - type(Mesh), intent(in) :: M - integer, intent(in) :: n - real(kind=rk), intent(in), target :: x(n) - real(kind=rk), intent(out), target :: r(n) - - real(kind=rk), pointer :: x_(:), r_(:) - - x_ => x - r_ => r - call SpMtx_pmvm(r_,A,x_,M) - -end subroutine ext_SpMtx_pmvm - -subroutine ext_vect_dot(v,x,y,n) - use Vect_mod - - real(kind=rk), intent(in), target :: x(n), y(n) - integer, intent(in) :: n - real(kind=rk), intent(out) :: v - - real(kind=rk), pointer :: x_(:), y_(:) - - x_ => x - y_ => y - v = Vect_dot_product(x_,y_) - -end subroutine ext_vect_dot - -subroutine ext_preconditioner_1level(sol,A,rhs,M,A_interf_,refactor_,nlf) - use pcg_mod - - implicit none - real(kind=rk),dimension(nlf),target :: sol !< solution - type(SpMtx) :: A !< sparse system matrix - real(kind=rk),dimension(nlf),target :: rhs !< right hand side - type(Mesh),intent(in) :: M !< Mesh - type(SpMtx),optional :: A_interf_ !< matr@interf. - logical,intent(inout) :: refactor_ - integer, intent(in) :: nlf - - !< residual vector, allocated - !! here for multiplicative Schwarz - real(kind=rk),dimension(:),pointer :: res => NULL() - - real(kind=rk), pointer :: sol_(:), rhs_(:) - - sol_ => sol - rhs_ => rhs - sol_ = 0 - call preconditioner_1level(sol_,A,rhs_,M,res,A_interf_,refactor_) - - if (sctls%method/=0) then - call Add_common_interf(sol_,A,M) - endif - -end subroutine ext_preconditioner_1level Index: c/ext/doug_ext.h =================================================================== --- src/ext/doug_ext.h (revision a10a43b3ff2bc776d5dd891961637125c8f5cec1) +++ (revision ) @@ -1,15 +1,0 @@ -#ifndef DOUG_EXT_H -#define DOUG_EXT_H - -extern void ext_doug_init_(int *init_type); -extern void ext_doug_finalize_(); -extern void alloc_spmtx_(void **A); -extern void alloc_mesh_(void **M); -extern void ext_paralleldistributeinput_(void *M, void *A, double *b, int *nparts, void *A_ghost, int *nlf); -extern void ext_cg_(void *A, double *b, void *M, double *xl, int *nlf); -extern void ext_spmtx_pmvm_(double *b, void *A, double *xl, void *M, int *nlf); -extern double ext_vec_dot_(double *v, double *x, double *y, int *nlf); -extern void ext_pmvmcommstructs_init_(void *A, void *M); -extern void ext_preconditioner_1level_(double *sol, void *A, double *rhs, void *M, void *A_interf_, int *refactor_, int *nlf); - -#endif Index: c/ext/main_aggr.c =================================================================== --- src/ext/main_aggr.c (revision cbd442c8070f6cedf638d0fd7e3de5cd373a0dfd) +++ (revision ) @@ -1,37 +1,0 @@ -#include -#include - -#include "doug_ext.h" - -int main(int argc, char **argv) -{ - int i; - void *A, *A_ghost; - void *M; - double *xl; - double *b; - int nparts=1; - int init_type=1; - int nlf; - - xl = malloc(225*sizeof(double)); - b = malloc(225*sizeof(double)); - - ext_doug_init_(&init_type); // parallel - alloc_spmtx_(&A); - alloc_spmtx_(&A_ghost); - alloc_mesh_(&M); - printf("Initialized, matrix %x, mesh %x\n", A, M); - ext_paralleldistributeinput_(M, A, b, &nparts, A_ghost, &nlf); - - ext_cg_(A, b, M, xl, &nlf); - - ext_doug_finalize_(); - - printf("result size=%d:\n", nlf); - for(i=0; i -#include -#include - -#include - -#include "doug_ext.h" -#include "ops.h" - -static void cg(void *A, void *M, double *xl, double *b, int nlf) -{ - double *r = malloc(nlf*sizeof(double)); - double *p = malloc(nlf*sizeof(double)); - double *q = malloc(nlf*sizeof(double)); - double tol, rho, rho_old, alpha, beta, tmp; - int it; - int rank, i; - - tol = 1E-12; - - // initialize pmvm constructs for communication - ext_pmvmcommstructs_init_(A, M); - - rho_old = 1.0; - ext_spmtx_pmvm_(r,A,xl,M, &nlf); - daxpy(r, -1.0, r, b, nlf); - for(i=0; i tol*tol) { - ext_vect_dot_(&rho, r, r, &nlf); - if (rank==0) - printf("rank=%d, it = %d, sqrt(rho) = %e\n", rank, it, sqrt(rho)); - - beta = rho/rho_old; - daxpy(p, beta, p, r, nlf); - ext_spmtx_pmvm_(q,A,p,M, &nlf); - ext_vect_dot_(&tmp, p,q,&nlf); - alpha = rho/tmp; - daxpy(xl, alpha, p, xl, nlf); - daxpy(r, -alpha, q, r, nlf); - rho_old = rho; - it++; - } -} - -int main(int argc, char **argv) -{ - int i; - void *A, *A_ghost; - void *M; - double *xl; - double *b; - int nparts=1; - int init_type=1; - int nlf; - - xl = malloc(225*sizeof(double)); - b = malloc(225*sizeof(double)); - - ext_doug_init_(&init_type); // parallel - alloc_spmtx_(&A); - alloc_spmtx_(&A_ghost); - alloc_mesh_(&M); - printf("Initialized, matrix %x, mesh %x\n", A, M); - ext_paralleldistributeinput_(M, A, b, &nparts, A_ghost, &nlf); - - cg(A, M, xl, b, nlf); - - ext_doug_finalize_(); - - for(i=0; i -#include -#include - -#include - -#include "doug_ext.h" -#include "ops.h" - -static void pcg(void *A, void *M, void *A_ghost, double *xl, double *b, int nlf) -{ - double *r = malloc(nlf*sizeof(double)); - double *p = malloc(nlf*sizeof(double)); - double *q = malloc(nlf*sizeof(double)); - double *z = malloc(nlf*sizeof(double)); - double tol, rho, rho_old, alpha, beta, tmp; - double ratio_norm, init_norm, res_norm; - int it, refactor=1; - int rank, i; - - tol = 1E-12; - - // initialize pmvm constructs for communication - ext_pmvmcommstructs_init_(A, M); - - rho_old = 1.0; - ext_spmtx_pmvm_(r,A,xl,M, &nlf); - daxpy(r, -1.0, r, b, nlf); - - MPI_Comm_rank(MPI_COMM_WORLD, &rank); - - ext_vect_dot_(&init_norm, r,r,&nlf); - if (init_norm == 0.0) init_norm = 1.0; - ratio_norm = 1.0; - - it = 0; - while (ratio_norm > tol*tol) { - it++; - ext_preconditioner_1level_(z,A,r,M,A_ghost,&refactor,&nlf); - refactor = 0; - - ext_vect_dot_(&rho, r, z, &nlf); - - if (it==1) { - ops_copy(p, z, nlf); - } else { - beta = rho/rho_old; - daxpy(p, beta, p, z, nlf); - } - ext_spmtx_pmvm_(q,A,p,M, &nlf); - ext_vect_dot_(&tmp, p,q,&nlf); - alpha = rho/tmp; - daxpy(xl, alpha, p, xl, nlf); - daxpy(r, -alpha, q, r, nlf); - rho_old = rho; - - ext_vect_dot_(&res_norm, r, r, &nlf); - ratio_norm = res_norm / init_norm; - if (rank==0) - printf("rank=%d, it = %d, sqrt(res_norm) = %e\n", rank, it, sqrt(res_norm)); - } -} - -int main(int argc, char **argv) -{ - int i; - void *A, *A_ghost; - void *M; - double *xl; - double *b; - int nparts=1; - int init_type=1; - int nlf; - - xl = malloc(225*sizeof(double)); - b = malloc(225*sizeof(double)); - - ext_doug_init_(&init_type); // parallel - alloc_spmtx_(&A); - alloc_spmtx_(&A_ghost); - alloc_mesh_(&M); - printf("Initialized, matrix %x, mesh %x\n", A, M); - ext_paralleldistributeinput_(M, A, b, &nparts, A_ghost, &nlf); - - pcg(A, M, A_ghost, xl, b, nlf); - - ext_doug_finalize_(); - - for(i=0; i main program that accepts mesh from 'triangle' output +program main_triangle + use DOUG_utils + use Mesh_class + use read_triangle_data_mod + use triangle_mod + implicit none + + character (len=FN_LEN) :: filename = 'mesh2d.1' + type(triangle_data) :: triData + type(Mesh) :: Msh !< Mesh + integer :: nparts !< number of partitons to partition a mesh + !> partition options (see METIS manual) + integer, dimension(6) :: part_opts = (/0,0,0,0,0,0/) + + ! Init DOUG + call DOUG_Init() + + nparts = numprocs + + ! Read triangle data and transform to Mesh data + if (ismaster()) then + call read_triangle_data(filename, triData) + call triangle_to_mesh(triData, Msh) + endif + + ! Get from master Mesh's parameters: nell, ngf, mfrelt, nsd, nnode + call Mesh_paramsMPIBCAST(Msh) + ! Master non-blockingly sends mesh data: nfrelt, mhead + call Mesh_dataMPIISENDRECV(Msh, & + nfrelt = .true., & + mhead = .true.) + + ! Build dual graph + call Mesh_buildGraphDual(Msh) + + ! Partition mesh's dual graph + if (ismaster()) then + if (sctls%plotting == D_PLOT_YES) then + call Mesh_pl2D_mesh(Msh) + end if + + ! Partition graph: D_PART_PMETIS, D_PART_KMETIS, D_PART_VKMETIS + call Mesh_partitionDual(Msh, nparts, D_PART_VKMETIS, part_opts) + if (sctls%plotting == D_PLOT_YES) then + call Mesh_pl2D_partitions(Msh) + end if + endif + + ! Distribute elements to partitons map among slaves + call Mesh_dataMPIISENDRECV(Msh, eptnmap=.true.) + + call DOUG_Finalize() +end program main_triangle Index: src/ext/triangle/read_triangle_data_mod.f90 =================================================================== --- src/ext/triangle/read_triangle_data_mod.f90 (revision 4bce0762b411f260032c6e4d1f6bf6f607e6ef29) +++ src/ext/triangle/read_triangle_data_mod.f90 (revision 4bce0762b411f260032c6e4d1f6bf6f607e6ef29) @@ -0,0 +1,647 @@ +module read_triangle_data_mod + + !----------------------------------------------------------------------------! + ! BRIEF DESCRIPTION OF THE MODULE + !----------------------------------------------------------------------------! + ! + ! Contains the subroutine READ_TRIANGLE_DATA which reads output files with + ! extensions .node, .ele, .neigh and .edge of the 2d grid generator TRIANGLE + ! into a fortran derived type TRIANGLE_DATA which can be subsequently used in + ! a fortran code. + ! + ! Example of usage: + ! ================= + ! + ! Assuming we have four triangle output files named + ! + ! filename.node + ! filename.ele + ! filename.neigh + ! filename.edge + ! + ! we can copy all the data in these files into the fortran derived type TD + ! using the following call + ! + ! ... + ! use read_triangle_data_mod + ! + ! type (triangle_data) :: td + ! + ! call read_triangle_data( 'filename', & ! in + ! td ) ! out + ! ... + ! + !----------------------------------------------------------------------------! + ! + ! Adapted with slight modifications from a subroutine of the public code + ! PHAML version 1.8 (see http://math.nist.gov/phaml/) by Artan Qerushi. + ! + ! Last modified on 11/17/2010 + !----------------------------------------------------------------------------! + + implicit none + + private + public :: rk, point, triangle_data, FN_LEN + public :: read_triangle_data + + ! kind parameter + integer, parameter :: rk = selected_real_kind(15, 307) + + type :: point + real (kind=rk) :: x, y + end type point + + type :: triangle_data + type (point), dimension(:), pointer :: vert_coord + real (kind=rk), dimension(:), pointer :: vert_bparam + integer, dimension(:,:), pointer :: tri_edge, tri_vert, tri_neigh, edge_tri, & + edge_vert, vert_tri, vert_edge + integer, dimension(:), pointer :: edge_bmark, vert_bmark + integer :: ntri, nedge, nvert + end type triangle_data + + ! Caesar version number + integer, parameter :: version_number = 1.0 + + ! user input error parameter + integer, parameter :: USER_INPUT_ERROR = -2 + + integer, parameter :: FN_LEN = 256 + + ! RESTRICTION no more than 16 triangles share a vertex in the triangle data + integer, parameter :: MAX_TD_VERT_NEIGH = 16 + + integer, parameter :: VERTICES_PER_ELEMENT = 3, & + EDGES_PER_ELEMENT = 3 + + integer :: ierr + +contains + + !============================================================================! + + subroutine read_triangle_data( triangle_files, & + td ) + + !--------------------------------------------------------------------------! + ! This routine reads data from .node, .ele, .edge and .neigh files in the + ! format of Jonathan Richard Shewchuk's mesh generation program "triangle". + ! + ! NOTE: I assume there are no comment lines before the end of the data. + ! Triangle 1.5 seems to obey this. + !--------------------------------------------------------------------------! + + ! Arguments + character (len=FN_LEN), intent(in) :: triangle_files + type (triangle_data), intent(out) :: td + + ! Local variables: + integer :: i, j, k, stat, td_dim, td_natt, td_nbm, bmark, vert, td_npt, tri, & + v1, v2, v3, td_ntri2, td_nneigh, edge, end1, end2, iounit, nset + logical, dimension(:), allocatable :: used + real (kind=rk) :: x, y + logical :: exists, opened + + !--------------------------------------------------------------------------! + ! Begin executable code + + ! find an available i/o unit number + + iounit = 11 + + do + + inquire( unit = iounit, & + exist = exists, & + opened = opened ) + + if ( exists .and. .not. opened ) exit + + iounit = iounit + 1 + + end do + + ! read the node (vertex) data from the triangle data file + + open( unit = iounit, & + file = trim(triangle_files)//".node", & + status = "old", & + action = "read", & + iostat = stat ) + + if ( stat /= 0 ) then + + call fatal( "open failed for file "//trim(triangle_files)//".node", & ! + "iostat is ", & ! + intlist = (/ stat /) ) ! + + stop + + endif + + read(iounit,*) td % nvert, td_dim, td_natt, td_nbm + + if ( td_natt /= 0 ) then + + call fatal("number of attributes in .node file must be 0") + + stop + + endif + + allocate( used(td % nvert), & + stat = stat ) + + if ( stat /= 0 ) then + + call fatal("memory allocation for vertices from .node file failed", & + intlist = (/ stat, td % nvert /) ) + + stop + endif + + allocate( td % vert_tri(MAX_TD_VERT_NEIGH,td % nvert), & + td % vert_edge(MAX_TD_VERT_NEIGH,td % nvert), & + td % vert_bmark(td % nvert), td % vert_bparam(td % nvert), & + td % vert_coord(td % nvert), & + stat = stat ) + + if ( stat /= 0 ) then + + call fatal("memory allocation for vertices from .node file failed", & + intlist = (/ stat, td % nvert /) ) + + stop + + endif + + if ( td_nbm == 0 ) then + + call fatal("boundary markers are required in data from Triangle") + + stop + + endif + + do i = 1, td % nvert + + read(iounit,*) vert,x,y,bmark + + if ( vert < 1 ) then + + call fatal("vertices in .node file must be numbered starting at 1") + + stop + + endif + + if ( vert > td % nvert ) then + + call fatal("vertex number in .node file is larger than stated number of vertices", & + intlist = (/ vert, td % nvert /) ) + + stop + + endif + + td % vert_coord(vert) % x = x + td % vert_coord(vert) % y = y + td % vert_bmark(vert) = bmark + + end do + + close(unit=iounit) + + ! read the element data from the .ele file + + open( unit = iounit, & + file = trim(triangle_files)//".ele", & + status = "old", & + action = "read", & + iostat = stat ) + + if ( stat /= 0 ) then + + call fatal("open failed for file "//trim(triangle_files)//".ele", & + "iostat is ", & + intlist = (/ stat /) ) + + stop + + endif + + read(iounit,*) td % ntri, td_npt, td_natt + + allocate( td % tri_edge(EDGES_PER_ELEMENT,td % ntri), & + td % tri_vert(VERTICES_PER_ELEMENT,td % ntri), & + stat = stat ) + + if ( stat /= 0 ) then + call fatal("memory allocation for triangles from .ele file failed", & + intlist = (/ stat, td % ntri /) ) + stop + endif + + used = .false. + + do i = 1, td % ntri + + read(iounit,*) tri,v1,v2,v3 + + if ( tri < 1 ) then + + call fatal("triangles in .ele file must be numbered starting at 1") + stop + + endif + + if ( tri > td % ntri ) then + + call fatal("triangle number in .ele file is larger than stated number of triangles", & + intlist = (/ tri, td % ntri /) ) + + stop + + endif + + td % tri_vert(1,tri) = v1 + td % tri_vert(2,tri) = v2 + td % tri_vert(3,tri) = v3 + used(v1) = .true. + used(v2) = .true. + used(v3) = .true. + + end do + + close(unit=iounit) + + if ( .not. all(used) ) then + + ierr = USER_INPUT_ERROR + + call fatal("There are unused nodes in the .node file.", & + "Use the -j flag when running triangle.") + + stop + + endif + + deallocate(used) + + ! read the neighbor data from the .neigh file + + open( unit = iounit, & + file = trim(triangle_files)//".neigh", & + status = "old", & + action = "read", & + iostat = stat ) + + if ( stat /= 0 ) then + call fatal( "open failed for file "//trim(triangle_files)//".neigh", & + "iostat is ", & + intlist = (/ stat /) ) + stop + endif + + read(iounit,*) td_ntri2, td_nneigh + + if ( td_ntri2 /= td % ntri ) then + call fatal("number of triangles in .neigh file is not the same as number in .ele file", & + intlist = (/ td_ntri2, td % ntri /) ) + stop + + endif + + allocate( td % tri_neigh(3,td % ntri), & + stat = stat ) + + if ( stat /= 0 ) then + call fatal( "memory allocation for neighbors from .neigh file failed", & + intlist = (/ stat, td % ntri /) ) + + stop + + endif + + do i = 1, td % ntri + + read(iounit,*) tri,v1,v2,v3 + + if ( tri < 1 ) then + call fatal("triangles in .neigh file must be numbered starting at 1") + stop + endif + + if ( tri > td % ntri ) then + call fatal("triangle number in .neigh file is larger than stated number of triangles", & + intlist = (/ tri, td % ntri /) ) + stop + endif + + td % tri_neigh(1,tri) = v1 + td % tri_neigh(2,tri) = v2 + td % tri_neigh(3,tri) = v3 + + end do + + close(unit=iounit) + + ! read the edge data from the triangle edge file + + open( unit = iounit, & + file = trim(triangle_files)//".edge", & + status = "old", & + action = "read", & + iostat = stat ) + + if ( stat /= 0 ) then + + call fatal("open failed for file "//trim(triangle_files)//".edge", & + "iostat is ", & + intlist = (/ stat /) ) + + stop + + endif + + read(iounit,*) td % nedge, td_nbm + + allocate( td % edge_tri(2,td % nedge), & + td % edge_vert(2,td % nedge), & + td % edge_bmark(td % nedge), & + stat = stat ) + + if ( stat /= 0 ) then + + call fatal("memory allocation for vertices from .edge file failed", & + intlist = (/ stat, td % nedge /) ) + + stop + + endif + + if ( td_nbm == 0 ) then + + call fatal("boundary markers are required in data from Triangle") + + stop + + endif + + do i = 1, td % nedge + + read(iounit,*) edge,end1,end2,bmark + + if ( edge < 1 ) then + + call fatal("edges in .edge file must be numbered starting at 1") + + stop + + endif + + if ( edge > td % nedge ) then + + call fatal("edge number in .edge file is larger than stated number of edges", & + intlist = (/ edge, td % nedge /) ) + + stop + endif + + td % edge_vert(1,edge) = end1 + td % edge_vert(2,edge) = end2 + td % edge_bmark(edge) = bmark + + end do + + close(unit=iounit) + + ! NEW + ! derive other components of triangle data + + ! set the triangle list for each vertex + + td % vert_tri = -1 + + do i = 1, td % ntri + do j = 1, 3 + + do k = 1, MAX_TD_VERT_NEIGH + if ( td % vert_tri(k,td % tri_vert(j,i)) == -1 ) exit + end do + + if ( k == MAX_TD_VERT_NEIGH + 1 ) then + call fatal("too many neighbors of a vertex in triangle data") + stop + endif + + td % vert_tri(k,td % tri_vert(j,i)) = i + + end do + + end do + + ! set the edge list for each vertex + + td % vert_edge = -1 + + do i = 1, td % nedge + do j = 1, 2 + + do k = 1, MAX_TD_VERT_NEIGH + if ( td % vert_edge(k,td % edge_vert(j,i)) == -1 ) exit + end do + + if ( k == MAX_TD_VERT_NEIGH + 1 ) then + call fatal("too many neighbors of a vertex in triangle data") + stop + endif + + td % vert_edge(k,td % edge_vert(j,i)) = i + + end do + + end do + + ! set the edge list of each triangle, and triangle list of each edge + + td % tri_edge = -1 + td % edge_tri = -1 + + ! for each triangle + + do i = 1, td % ntri + + nset = 0 + + ! for each vertex of the triangle + + do j = 1, 3 + + ! search the edges of the vertex for any that contain another vertex of the + ! triangle + ! for each edge of this vertex + + do k = 1, MAX_TD_VERT_NEIGH + + if (td % vert_edge(k,td % tri_vert(j,i)) == -1) exit + + ! if the first vertex of the edge is this vertex, see if the second vertex of + ! the edge is a vertex of the triangle + + if ( td % edge_vert(1,td % vert_edge(k,td % tri_vert(j,i))) == td % tri_vert(j,i) ) then + + if ( td % edge_vert(2,td % vert_edge(k,td % tri_vert(j,i)) ) == & + td % tri_vert(1,i) .or. & + td % edge_vert(2,td % vert_edge(k,td % tri_vert(j,i))) == & + td % tri_vert(2,i) .or. & + td % edge_vert(2,td % vert_edge(k,td % tri_vert(j,i))) == & + td % tri_vert(3,i) ) then + + ! if so make it an edge of this triangle, and make this triangle a triangle + ! of that edge, unless it has already been set + + if ( td % edge_tri(1,td % vert_edge(k,td % tri_vert(j,i)) ) /= i .and. & + td % edge_tri(2,td % vert_edge(k,td % tri_vert(j,i)) ) /= i ) then + + nset = nset + 1 + td % tri_edge(nset,i) = td % vert_edge(k,td % tri_vert(j,i)) + + if ( td % edge_tri(1,td % vert_edge(k,td % tri_vert(j,i))) == -1 ) then + + td % edge_tri(1,td % vert_edge(k,td % tri_vert(j,i))) = i + + elseif ( td % edge_tri(2,td % vert_edge(k,td % tri_vert(j,i))) == -1 ) then + + td % edge_tri(2,td % vert_edge(k,td % tri_vert(j,i))) = i + + else + + call fatal("too many triangles neighboring an edge in read_trianlge_data") + + stop + + endif + + endif + + endif + + ! if the second vertex of the edge is this vertex, see if the first vertex of + ! the edge is a vertex of the triangle + + elseif ( td % edge_vert(2,td % vert_edge(k,td % tri_vert(j,i))) == td % tri_vert(j,i) ) then + + if ( td % edge_vert(1,td % vert_edge(k,td % tri_vert(j,i))) == & + td % tri_vert(1,i) .or. & + td % edge_vert(1,td % vert_edge(k,td % tri_vert(j,i))) == & + td % tri_vert(2,i) .or. & + td % edge_vert(1,td % vert_edge(k,td % tri_vert(j,i))) == & + td % tri_vert(3,i) ) then + + ! if so make it an edge of this triangle, and make this triangle a triangle + ! of that edge, unless it has already been set + + if ( td % edge_tri(1,td % vert_edge(k,td % tri_vert(j,i))) /= i .and. & + td % edge_tri(2,td % vert_edge(k,td % tri_vert(j,i))) /= i ) then + + nset = nset + 1 + td % tri_edge(nset,i) = td % vert_edge(k,td % tri_vert(j,i)) + + if (td % edge_tri(1,td % vert_edge(k,td % tri_vert(j,i))) == -1) then + + td % edge_tri(1,td % vert_edge(k,td % tri_vert(j,i))) = i + + elseif (td % edge_tri(2,td % vert_edge(k,td % tri_vert(j,i))) == -1) then + + td % edge_tri(2,td % vert_edge(k,td % tri_vert(j,i))) = i + + else + + call fatal("too many triangles neighboring an edge in read_trianlge_data") + + stop + + endif + + endif + + endif + + endif + + end do + + end do + + end do + + ! verify that all triangles have 3 edges and all edges have 2 triangles or + ! are boundary + + do i = 1, td % ntri + if ( td % tri_edge(3,i) == -1 ) then + + call fatal("didn't assign 3 edges to all triangles in read_triangle_data") + + stop + + endif + end do + + ! TEMP must verify that bmark/=0 iff vertex is on boundary. Might need to + ! change the documentation + + do i = 1, td % nedge + if ( td % edge_tri(1,i) == -1 .or. & + ( td % edge_bmark(i) == 0 .and. td % edge_tri(2,i) == -1 ) ) then + + call fatal("didn't assign 2 triangles or 1 triangle and boundary mark to all edges in read_triangle_data") + + stop + + endif + end do + + end subroutine read_triangle_data + + !============================================================================! + + subroutine fatal( msg, & + msg2, & + intlist, & + reallist ) + + !--------------------------------------------------------------------------! + ! This routine handles fatal errors + !--------------------------------------------------------------------------! + + ! Arguments + character (len=*), intent(in) :: msg + character (len=*), intent(in), optional :: msg2 + integer, intent(in), optional, dimension(:) :: intlist + real (kind=rk), intent(in), optional, dimension(:) :: reallist + + !--------------------------------------------------------------------------! + ! Begin executable code + + write(*,"(A)") + write(*,"(A)") "------------------------------------------------------" + write(*,"(3A)") " Caesar Version ", version_number, " ERROR" + write(*,"(A)") msg + + if ( present(msg2) ) write(*,"(A)") msg2 + if ( present(intlist) ) write(*,"(7I11)") intlist + if ( present(reallist) ) write(*,"(SS,1P,4E18.10E2)") reallist + + write(*,"(A)") "------------------------------------------------------" + write(*,"(A)") + + stop + + end subroutine fatal + + !============================================================================! + +end module read_triangle_data_mod Index: src/ext/triangle/triangle_mod.f90 =================================================================== --- src/ext/triangle/triangle_mod.f90 (revision 4bce0762b411f260032c6e4d1f6bf6f607e6ef29) +++ src/ext/triangle/triangle_mod.f90 (revision 4bce0762b411f260032c6e4d1f6bf6f607e6ef29) @@ -0,0 +1,64 @@ +!> Module that handles data from 'triangle' program output +module triangle_mod + use Mesh_class + use read_triangle_data_mod + +contains + subroutine triangle_to_mesh(td, M) + implicit none + type(triangle_data), intent(in) :: td + type(Mesh), intent(out) :: M + + integer :: i,j,k + integer, allocatable :: all2freeMap(:) ! index map from all nodes to free nodes + + M = Mesh_New() + M%nell = td%ntri + M%nnode = td%nvert + M%ngf = count(td%vert_bmark==0)! number of free nodes (i.e. non-fixed nodes) <=nnode + M%mfrelt = 3 ! max number of free nodes per element (i.e. VERTICES_PER_ELEMENT) + M%nsd = 2 ! number of spatial dimensions + + ! create index map from free nodes to all nodes + allocate(M%freemap(M%ngf)) + k = 0 + do i = 1, M%nnode + if (td%vert_bmark(i)==0) then + k = k+1 + M%freemap(k) = i + end if + end do + + ! create reverse map + allocate(all2freeMap(M%nnode)) + all2freeMap = 0 + all2freeMap(M%freemap) = (/ (i,i=1,M%ngf) /) + + ! copy element info + allocate(M%nfrelt(M%nell)) ! number of free nodes for each element + allocate(M%mhead(M%mfrelt,M%nell)) ! free node indices for each element + + ! copy triangle vertices (DOUG has 0 as missing) + M%mhead = 0 + do i = 1, M%nell + ! calculate number of free (non-boundary) nodes + M%nfrelt(i) = count(td%vert_bmark(td%tri_vert(:,i))==0) + k = 0 + do j = 1, 3 ! VERTICES_PER_ELEMENT + ! exclude boundary nodes (treat as fixed) + if (td%vert_bmark(td%tri_vert(j,i)) == 0) then + k = k+1 + M%mhead(k,i) = all2freeMap(td%tri_vert(j,i)) + end if + end do + end do + + ! copy coordinates + allocate(M%coords(M%nsd,M%nnode)) + do i = 1, M%nnode + M%coords(1,i) = td%vert_coord(i)%x + M%coords(2,i) = td%vert_coord(i)%y + end do + + end subroutine triangle_to_mesh +end module triangle_mod Index: src/main/main_drivers.F90 =================================================================== --- src/main/main_drivers.F90 (revision 848affb82d75093fa2bc5d79fa514342ac979b36) +++ src/main/main_drivers.F90 (revision c28debaaa859fac8b6880876f9b655be422a21b5) @@ -88,4 +88,16 @@ type(SpMtx),intent(in out),optional :: A_interf !< matrix at interface + if (ismaster()) then ! MASTER + write(stream,*) + write(stream,*) 'master thread' + if (D_MSGLVL > 1) & + call MasterCtrlData_print() + + else ! SLAVES + write(stream,'(a,i4,a)') 'slave [',myrank,'] thread' + if (D_MSGLVL > 1) & + call SharedCtrlData_print() + end if + ! ======================= ! Mesh and its Dual Graph @@ -93,22 +105,11 @@ ! ! Create Mesh object + Msh = Mesh_New() - if (ismaster()) then ! MASTER - write(stream,*) - write(stream,*) 'master thread' - - if (D_MSGLVL > 1) & - call MasterCtrlData_print() - + if (ismaster()) then ! Initialise Mesh object call Mesh_initFromFile(Msh, trim(mctls%info_file)) - - else ! SLAVES - write(stream,'(a,i4,a)') 'slave [',myrank,'] thread' - - if (D_MSGLVL > 1) & - call SharedCtrlData_print() - end if + endif ! Get from master Mesh's parameters: nell, ngf, mfrelt, nsd, nnode @@ -116,11 +117,4 @@ if (D_MSGLVL > 1) & call Mesh_printInfo(Msh) - - ! Allocate data arrays (nfrelt, mhead, freemap, eptnmap) for mesh - call Mesh_allocate(Msh, & - nfrelt =.true., & - mhead =.true., & - freemap =.true., & - eptnmap =.true.) ! Master reads in from files: feedom lists, coordinates, freedom map @@ -141,5 +135,4 @@ ! For multi-variable problems which have more than one block if (sctls%number_of_blocks > 1) then - call Mesh_allocate(Msh, freemask=.true.) if (ismaster()) then call Mesh_readFromFile(Msh, & @@ -156,43 +149,14 @@ ! Partition mesh's dual graph if (ismaster()) then - - !! Build dual graph (Graph object is a data field in Mesh class) - !call Mesh_buildGraphDual(Msh) - if (sctls%plotting == D_PLOT_YES) then - - call Mesh_pl2D_pointCloud(Msh,D_PLPLOT_INIT) - ! Plots mesh's dual graph - call Mesh_pl2D_plotGraphDual(Msh,D_PLPLOT_END) - ! Mesh & its Dual Graph - call Mesh_pl2D_plotMesh(Msh, D_PLPLOT_INIT) - call Mesh_pl2D_plotGraphDual(Msh, D_PLPLOT_END) + call Mesh_pl2D_mesh(Msh) end if ! Partition graph: D_PART_PMETIS, D_PART_KMETIS, D_PART_VKMETIS call Mesh_partitionDual(Msh, nparts, D_PART_VKMETIS, part_opts) - if (sctls%plotting == D_PLOT_YES) then - ! Draw colored partitoined graph - call Mesh_pl2D_plotGraphParted(Msh) - - ! Plot partitions of the mesh - ! NB: Check for multivariable case! TODO - call Mesh_pl2D_Partition(Msh) - ! Partition with Dual Graph upon it - call Mesh_pl2D_Partition(Msh, D_PLPLOT_INIT) - call Mesh_pl2D_plotGraphDual(Msh, D_PLPLOT_CONT) - call Mesh_pl2D_pointCloud(Msh,D_PLPLOT_END) - end if - - ! Destroy previously created graph (purely to save memory) - ! (If it is not killed here or somewere else Mesh_Destroy() - ! will kill it any way) - !call Mesh_destroyGraph(Msh) - else ! SLAVES - ! Number of partions the mesh was partitioned into - Msh%nparts = nparts - Msh%parted = .true. - end if + call Mesh_pl2D_partitions(Msh) + end if + endif ! Distribute elements to partitons map among slaves