forums.silverfrost.com Forum Index forums.silverfrost.com
Welcome to the Silverfrost forums
 
 FAQFAQ   SearchSearch   MemberlistMemberlist   UsergroupsUsergroups   RegisterRegister 
 ProfileProfile   Log in to check your private messagesLog in to check your private messages   Log inLog in 

Simple 2D mesher (in Fortran)
Goto page Previous  1, 2
 
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> General
View previous topic :: View next topic  
Author Message
jjgermis



Joined: 21 Jun 2006
Posts: 404
Location: Nürnberg, Germany

PostPosted: Mon Aug 22, 2011 9:36 am    Post subject: Reply with quote

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 Very Happy 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
View user's profile Send private message
jjgermis



Joined: 21 Jun 2006
Posts: 404
Location: Nürnberg, Germany

PostPosted: Wed Aug 24, 2011 7:56 am    Post subject: Reply with quote

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
View user's profile Send private message
jjgermis



Joined: 21 Jun 2006
Posts: 404
Location: Nürnberg, Germany

PostPosted: Wed Aug 24, 2011 11:45 am    Post subject: Reply with quote

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
View user's profile Send private message
jjgermis



Joined: 21 Jun 2006
Posts: 404
Location: Nürnberg, Germany

PostPosted: Mon Aug 29, 2011 1:20 pm    Post subject: Reply with quote

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
View user's profile Send private message
DanRRight



Joined: 10 Mar 2008
Posts: 2139
Location: South Pole, Antarctica

PostPosted: Sun Sep 04, 2011 3:03 pm    Post subject: Reply with quote

Cool. Keep us updating on the progress. Great would be if you write some simple program examples for everyone. Any plans to move to 3D like in this link below?

https://computation.llnl.gov/casc/Overture/henshaw/publications/OvertureShortCourse.pdf
Back to top
View user's profile Send private message
jjgermis



Joined: 21 Jun 2006
Posts: 404
Location: Nürnberg, Germany

PostPosted: Tue Sep 06, 2011 12:20 pm    Post subject: Reply with quote

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 Smile

Back to top
View user's profile Send private message
jjgermis



Joined: 21 Jun 2006
Posts: 404
Location: Nürnberg, Germany

PostPosted: Wed Sep 07, 2011 12:25 pm    Post subject: Reply with quote

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
View user's profile Send private message
jjgermis



Joined: 21 Jun 2006
Posts: 404
Location: Nürnberg, Germany

PostPosted: Thu Sep 08, 2011 10:46 am    Post subject: Reply with quote

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
View user's profile Send private message
jjgermis



Joined: 21 Jun 2006
Posts: 404
Location: Nürnberg, Germany

PostPosted: Fri Sep 09, 2011 11:54 am    Post subject: Reply with quote

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
View user's profile Send private message
DanRRight



Joined: 10 Mar 2008
Posts: 2139
Location: South Pole, Antarctica

PostPosted: Sat Sep 10, 2011 10:11 pm    Post subject: Reply with quote

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
View user's profile Send private message
jjgermis



Joined: 21 Jun 2006
Posts: 404
Location: Nürnberg, Germany

PostPosted: Mon Sep 12, 2011 2:58 pm    Post subject: Reply with quote

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 Very Happy
Back to top
View user's profile Send private message
jjgermis



Joined: 21 Jun 2006
Posts: 404
Location: Nürnberg, Germany

PostPosted: Fri Sep 30, 2011 10:05 am    Post subject: Reply with quote

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
View user's profile Send private message
jjgermis



Joined: 21 Jun 2006
Posts: 404
Location: Nürnberg, Germany

PostPosted: Mon Mar 26, 2012 3:06 pm    Post subject: Reply with quote

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
View user's profile Send private message
Display posts from previous:   
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> General All times are GMT + 1 Hour
Goto page Previous  1, 2
Page 2 of 2

 
Jump to:  
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