|
forums.silverfrost.com Welcome to the Silverfrost forums
|
View previous topic :: View next topic |
Author |
Message |
jjgermis
Joined: 21 Jun 2006 Posts: 404 Location: Nürnberg, Germany
|
Posted: Mon Aug 22, 2011 9:36 am Post subject: |
|
|
The GEOMPACK (Fortran version) software package compiles without any problems using FTN95. Furthermore, the following documents give a good description of the package:
1.) Barry Joe, "GEOMPACK - a software package for the generation of meshes using geometric algorithms", Advances in Engineering Software, Vol. 13, Nr. 5/6, 1991.
2.) File format: The file formats are online available (including examples). The successor for Fortran version is GEOMPACK++ in C++. The file formats are the same.
The author claims that the C++ version is 1.5 to 2 times faster is than the Fortran version (mainly bacause of using classes). I always thouht that Fortran is famous for number crunching - and is mesh generation is nothing else than number crunching.
In a previous thread I got very useful tipps on using Clearwin+ and drawing polygons One of the provided examples is shown below. Most of the meshing programs I have tested use the basically the same data structure to present the data. This means I can keep my original method of calculating the basic coordinates and use the subroutines from GEOMPACK to generate the mesh.
Uploaded with ImageShack.us |
|
Back to top |
|
|
jjgermis
Joined: 21 Jun 2006 Posts: 404 Location: Nürnberg, Germany
|
Posted: Wed Aug 24, 2011 7:56 am Post subject: |
|
|
I fortgot to mention that there is one thing to change when compiling with FTN: change the kind values! The following (simple) code reads the region and mesh files. For simple polygons one can ignore the curve file.
Code: | module GEOMType
type rg2_type
integer :: nvc
double precision,dimension(1000) :: x,y
integer,dimension(1000) :: vertinfo
integer :: nvx
integer, dimension(1000) :: nodecode,icurv
double precision,dimension(1000) :: ucurv
integer :: nloop
integer, dimension(1000) :: nv,regcode,inreg
integer, dimension(100,100) :: v,edginfo
end type rg2_type
type mh2_type
integer :: nvc
double precision,dimension(1000) :: x,y
integer,dimension(1000) :: vertinfo
integer :: nvx
integer, dimension(1000) :: nodecode,icurv
double precision,dimension(1000) :: ucurv
integer :: nodelem, nelem
integer,dimension(4,1000) :: v
integer,dimension(1000) :: regcode
integer,dimension(4,1000) :: edginfo
end type mh2_type
contains
subroutine read_rg2(rg2_file,rg2)
implicit none
character(len=*), intent(in) :: rg2_file
type(rg2_type), intent(out) :: rg2
integer, parameter :: iu = 23
integer :: i,j
write ( *, '(a)' ) ' Read '//trim(rg2_file)
open(iu,file=rg2_file)
read(iu,*) rg2%nvc
read(iu,*) (rg2%x(i),rg2%y(i),rg2%vertinfo(i),i=1,rg2%nvc)
read(iu,*) rg2%nvx
read(iu,*) (rg2%nodecode(i),rg2%icurv(i),rg2%ucurv(i),i=1,rg2%nvx)
read(iu,*) rg2%nloop
do i=1,rg2%nloop
read(iu,*) rg2%nv(i),rg2%regcode(i),rg2%inreg(i)
read(iu,*) (rg2%v(j,i),j=1,rg2%nv(i))
read(iu,*) (rg2%edginfo(j,i),j=1,rg2%nv(i))
enddo
return
end subroutine read_rg2
subroutine read_mh2(mh2_file,mh2)
implicit none
character(len=*), intent(in) :: mh2_file
type(mh2_type), intent(out) :: mh2
integer, parameter :: iu = 24
integer :: i,j
write ( *, '(a)' ) ' Read '//trim(mh2_file)
open(iu,file=mh2_file)
read(iu,*) mh2%nvc
read(iu,*) (mh2%x(i),mh2%y(i),mh2%vertinfo(i),i=1,mh2%nvc)
read(iu,*) mh2%nvx
read(iu,*) (mh2%nodecode(i),mh2%icurv(i),mh2%ucurv(i),i=1,mh2%nvx)
read(iu,*) mh2%nodelem,mh2%nelem
read(iu,*) ((mh2%v(i,j),i=1,mh2%nodelem),j=1,mh2%nelem)
read(iu,*) (mh2%regcode(i),(mh2%edginfo(j,i),j=1,mh2%nodelem),i=1,mh2%nelem)
return
end subroutine read_mh2
end module GEOMType
|
|
|
Back to top |
|
|
jjgermis
Joined: 21 Jun 2006 Posts: 404 Location: Nürnberg, Germany
|
Posted: Wed Aug 24, 2011 11:45 am Post subject: |
|
|
Another interessting (Fortran) mesher which is referenced in [1] is explained in a paper by Sloan [2]. In the paper the code is listed as well - after digging it up in our library I found an online copy. Moreover, more details on the examples given in [1] can be found in [3].
1.) Barry Joe, "GEOMPACK - a software package for the generation of meshes using geometric algorithms", Advances in Engineering Software, Vol. 13, Nr. 5/6, 1991.
2.) S.W. Sloan, "A fast algorithm for constructing Delaunay triangulations in the plane", Advances in Engineering Software, Vol. 9, Nr. 1, 1987.
3.) Barry Joe, "Triangular meshes for regions of complicated shape", Int. journal for numerical methods in eng., Vol. 23, pp. 751-778, 1986.
This link has a postscript copy of [1].
For obtaining triangular meshes the examle test17 that comes with the package should be studied. The result is given in [1]. |
|
Back to top |
|
|
jjgermis
Joined: 21 Jun 2006 Posts: 404 Location: Nürnberg, Germany
|
Posted: Mon Aug 29, 2011 1:20 pm Post subject: |
|
|
I would like to conclude this thread with an example form the above mentioned papers.
The choice of mesher depends on the type of application and the effort one would like to invest in mesh generation. Our application requires automatic mesh generation, a minimum input paramters and it must be in Fortran. The mesher presented here:
1.) requiers the vertex coordinates; and
2.) mesh control parameters.
Furthermore, the examples provided in the papers give all the information that is necessary to reproduce the figures and results presented in the testcode. The regions that the package requires is much the same as we use before meshing, i.e. adapting our code is very easy.
|
|
Back to top |
|
|
DanRRight
Joined: 10 Mar 2008 Posts: 2816 Location: South Pole, Antarctica
|
|
Back to top |
|
|
jjgermis
Joined: 21 Jun 2006 Posts: 404 Location: Nürnberg, Germany
|
Posted: Tue Sep 06, 2011 12:20 pm Post subject: |
|
|
Hi Dan
Thanks for the positive response! The (LaTeX) slides with the examples looks very impressive.
The mesher is indeed very cool and fulfill my requirements, i.e. high quality mesh with little input. For the moment I only focus on 2D. However, the mesher includes 3D routines as well.
In the meshing business it seems like Lake Superior is a popular goal - Triangle uses the same example. Below is the result of the mesher under discussion (I still have to include the holes).
As you can see in the previous thread my focus is on rotating objects. For this we use a mesh free air gap. The two parts are then "linked" by applying the LaPlace equation to the uniform part of the air gap - this is only a few extra lines in the direct solver (once the Fourier coefficients are calculated for the uniform part).
I would like to share some examples. However, in this regard this forum has potential
|
|
Back to top |
|
|
jjgermis
Joined: 21 Jun 2006 Posts: 404 Location: Nürnberg, Germany
|
Posted: Wed Sep 07, 2011 12:25 pm Post subject: |
|
|
Even when including the holes the mesher (almost) gets the job done. For some reason I get some problems with the tolerances (code below). Assume it must have somthing to do the subroutine where the angles are calculated - always a tricky part. However, when I remove the test, then the result is fine.
Code: | iang(vr) = angle ( xlft, ylft, xt, yt, xrgt, yrgt )
if ( iang(vr) < pi - tol .or. pi + tol < iang(vr) ) then
ierror = 219
return
end if |
Uploaded with ImageShack.us |
|
Back to top |
|
|
jjgermis
Joined: 21 Jun 2006 Posts: 404 Location: Nürnberg, Germany
|
Posted: Thu Sep 08, 2011 10:46 am Post subject: |
|
|
The problem with the tolerance was solved by introducing a global variable tol instead of local variables. This can be set in the input file. Calculating the tolerance is done as follows:
Code: | subroutine initcb(tolin)
implicit none
double precision,intent(out) :: tolin
double precision :: eps,epsp1
eps = 1.0d0
epsp1 = 1.5d0
do while (epsp1 > 1.0d0)
eps = eps/2.0d0
epsp1 = 1.0d0+eps
enddo
tol = max(tolin,100d0*eps)
pi = 4.0d0*datan(1.0d0)
return
end subroutine initcb |
The advantage of this software package is that the number of elements can be controlled by a single variable. This makes it very user friendly. Comparing with the initial mesh (see thread at the begining) one now can control the mesh. The original sotware is unchanged!
Uploaded with ImageShack.us |
|
Back to top |
|
|
jjgermis
Joined: 21 Jun 2006 Posts: 404 Location: Nürnberg, Germany
|
Posted: Fri Sep 09, 2011 11:54 am Post subject: |
|
|
Working through the papers by Barry Joe one should be familiar with the following definitions:
1.) simple polygon
2.) multiply connected region
3.) convex polygon
The principle of decomposition on which the GEOMPACK software package is based makes it very easy to automatically detect the boundaries.
Uploaded with ImageShack.us |
|
Back to top |
|
|
DanRRight
Joined: 10 Mar 2008 Posts: 2816 Location: South Pole, Antarctica
|
Posted: Sat Sep 10, 2011 10:11 pm Post subject: |
|
|
Well...Still you clearly see both mesh generators are creating not the meshes you'd call perfect. Both are very inhomogeneous. First one also creates very elongated triangles, the second one fixes that but the local inhomogeneity of density of triangles is still clearly visible. There exist not much of physical reason for that high inhomogeneity.
Is there any "random generator" mesh generators which place the points quasi homogeneously and then this mesh is improved based on the similar principles like GEOMPACK? I suspect they'd do the job much better in this at least respect. They could be slower though |
|
Back to top |
|
|
jjgermis
Joined: 21 Jun 2006 Posts: 404 Location: Nürnberg, Germany
|
Posted: Mon Sep 12, 2011 2:58 pm Post subject: |
|
|
All your comments on the inhomegeneous meshes are correct. There are several good mesh generators available. However, the only good ones in public domain (or known to me) are in C/C++. GEOMPACK ist written in Fortran.
On the other hand one should distinguish between meshes for mechanical and electrical applications. We have experienced in the past that a mesh for a mechanical problem is impractical for an electrical problem. The calculated flux linkages for the above meshes are pretty much the same.
Furthermore, we wish to integrate FEM in an everyday design environment, i.e. for a given machine topology the whole process must be fully automated. This can only be achieved if the mesher does a good job. For this purpose GEOMPACK does a fairly good job. The user only have to specify the number of triangles that is required and a smoothing factor. The latter controls the changeover between the "decomposed" simple convex polygons. Moreover, GEOMPACK automatically devide the periodic boundaries in such a manner that inserted nodes on the periodic boundaries maps.
All in all one have to compromise between what is nice to have and what is available. The for the rotating objects (and mesh free air gap) is written in Fortran. Even though mix programming is no problem I prefer a homogenous source code |
|
Back to top |
|
|
jjgermis
Joined: 21 Jun 2006 Posts: 404 Location: Nürnberg, Germany
|
Posted: Fri Sep 30, 2011 10:05 am Post subject: |
|
|
Using GEOMPACK as basis it is relative easy to adapt it for a specific application. A very good property is that the algorithm automatically determines the length scale in a specific region. This assures a good mesh quality overall and at cross overs between regions.
This mesher requires a definition of the regions which define the modell. For this the basic coordinates needs to be defined. Adding nodes afterwards means the region (variables) needs to be redefined. Due to the mesh variables this is no problem keeping track of the added nodes. The subroutine below adds nodes and update the mesh variables autmatically.
John: In the meanwhile I have learned a lot about meshing and produced quite a few lines of code. The whole time I was keeping the advancing front method in mind. Your hint motivated me to keep this in mind and I think that a some stage I will try to implement it!
Code: | subroutine nc_line(a,b,segnr,M,srg)
implicit none
integer, intent(in) :: a,b,segnr,M
type(srg_type) :: srg
integer :: i,j,n,ib,ie,d,e,k,ivrt_max
double precision :: dx,dy
!
! Add n-1 nodes to the segment with nodes from numbered d to e. First:
! 1.) find the starting node in the sequence;
! 2.) shift the nodes from k+1 to ivrt_max by n-1 positions to the right
! 3.) add the new nodes to the list
! 4.) update srg with new added nodes
!
dx = srg%x(a)-srg%x(b)
dy = srg%y(a)-srg%y(b)
n = int( dsqrt(dx*dx+dy*dy)/(srg%slot_pitch/dble(M)) + 0.5d0)
call calculate_nvc(srg)
ivrt_max = sum(srg%nvbc(1:srg%ncur))
d = sum(srg%nvbc(1:segnr-1))+1
e = d+srg%nvbc(segnr)-1
do j=d,e
if ( iabs(srg%ivrt(j)) == a ) then
k = j
exit
endif
enddo
ib = k+1+n-1
ie = ivrt_max+n-1
srg%ivrt(ib:ie) = srg%ivrt(k+1:ivrt_max)
srg%ivrt(k+1:e) = 0
dx = (srg%x(b)-srg%x(a))/dble(n)
dy = (srg%y(b)-srg%y(a))/dble(n)
do i = k+1,ib-1
srg%nvc = srg%nvc+1
srg%x(srg%nvc) = srg%x(a)+dble(i-k)*dx
srg%y(srg%nvc) = srg%y(a)+dble(i-k)*dy
srg%ivrt(i) = sign(srg%nvc,srg%ivrt(k))
end do
srg%nvc = srg%nvc+(n-1)
srg%nvbc(segnr) = srg%nvbc(segnr)+(n-1)
return
end subroutine |
|
|
Back to top |
|
|
jjgermis
Joined: 21 Jun 2006 Posts: 404 Location: Nürnberg, Germany
|
Posted: Mon Mar 26, 2012 3:06 pm Post subject: |
|
|
As already mentioned: The Goempack mesher works very well and offers good examples to test and for getting started. However, the documentation presents the required info very scattered. Working through the docs the following testvectors covers more than enough to fully make use of the mesh functions:
|
|
Back to top |
|
|
|
|
You cannot post new topics in this forum You cannot reply to topics in this forum You cannot edit your posts in this forum You cannot delete your posts in this forum You cannot vote in polls in this forum
|
Powered by phpBB © 2001, 2005 phpBB Group
|