Simple 2D mesher (in Fortran)

We use a very simple 2D mesh generator. The basic code is given below. It only takes a polygon, calculate the geometric mean and connect it with the nodes on the edges - very simple. A very good mesh can be found at triangle which compiles with SCC. Searching for meshes results in a never ending story... I would appreciate some tipps for a 2D mesher in Fortran.

    ! READ THE FILE CONTAINING THE POLYGON INFORMATION
    pol_file = TRIM(fls%dir_pol)//fls%SEP//TRIM(fls%pol_file)
    write(*,*) 'Read: '//TRIM(pol_file)
    open(unit=iu,FILE=pol_file)
    rewind iu
    read (iu,*,err=81) npol
    do i=1,npol
        read (iu,*,err=81) nspnts_pol,nstype_pol
        xint = 0.0d0
        yint = 0.0d0
        do j=1,nspnts_pol
            read (iu,*,err=81) x_ply(j),y_ply(j)
            xint = xint+x_ply(j)
            yint = yint+y_ply(j)
        end do 
    !   Calculate internal point of the polygon segment
        xint = xint/dble(nspnts_pol)
        yint = yint/dble(nspnts_pol)
        if (fpl%nnode == 0) then
        !     Assign the first node
            fpl%nnode    = 1
            fpl%x(1)     = x_ply(1)
            fpl%y(1)     = y_ply(1)
            ply_pnter(1) = 1
        end if
        do j=1,nspnts_pol
            do k=1,fpl%nnode
                dx   = x_ply(j)-fpl%x(k)
                dy   = y_ply(j)-fpl%y(k)
                dist = (dx*dx+dy*dy)**0.5d0
                if (dist.lt.told) then
                    ply_pnter(j) = k
                    goto 43
                end if
            end do
            fpl%nnode    = fpl%nnode+1
            ply_pnter(j) = fpl%nnode
            fpl%x(fpl%nnode) = x_ply(j)
            fpl%y(fpl%nnode) = y_ply(j)
            43 continue
        end do
        do k=1,fpl%nnode
            dx   = dabs(xint-fpl%x(k))
            dy   = dabs(yint-fpl%y(k))
            dist = (dx*dx+dy*dy)**0.5d0
            if (dist < told) then
                int_pnter(i) = k
                goto 44
            end if
        end do
        fpl%nnode    = fpl%nnode+1
        int_pnter(i) = fpl%nnode
        fpl%x(fpl%nnode) = xint
        fpl%y(fpl%nnode) = yint
        44 continue
        do j=1,(nspnts_pol-1)
            fpl%nelm     = fpl%nelm+1
            fpl%node(fpl%nelm,1) = ply_pnter(j)
            fpl%node(fpl%nelm,2) = ply_pnter((j+1))
            fpl%node(fpl%nelm,3) = int_pnter(i)
            fpl%itype(fpl%nelm)  = nstype_pol
        end do
        fpl%nelm         = fpl%nelm+1
        fpl%node(fpl%nelm,1) = ply_pnter(nspnts_pol)
        fpl%node(fpl%nelm,2) = ply_pnter(1)
        fpl%node(fpl%nelm,3) = int_pnter(i)
        fpl%itype(fpl%nelm)  = nstype_pol
    end do
    close(unit=iu,status='keep')

Hello Jjgermis,

I suggest that you try writing your own 2D mesher using the advancing front method. The logic of this is easy to understand and implement. When this is complete use a Laplacian smoothing algorithm to improve the element quality.

Good luck !

Hi John,

thanks for the tipp!

The key words advancing front method is already a good starting point. As with my resent literature study on direct solvers and renumbering I made good progress once I got the correct terminology.

Hi John, attached is a result of our current mesher. My experience for flux linkage calculations is that a coarse mesh is more than efficient and fast. Flux densities on the othet hand requires a much finer mesh. As you can see we use no mesh in the air gap.

FEM literature seems to be very limited for electrical applications. A good despription can be found in the book by Silverter and Ferrari. However, they only provide the theory and some (unoptimised) example code. For real software implementations I found the solvers and efficient use of memory in mechanical engineering FEM!

I think that FEM software is as good as its mesher! The solving part is relative simple. Obtaining a (quality) mesh is the name of the game!!

Following your tipp on advancing front I got a good idea and a few interessting websites. I think it might be possible to keep our existing code and apply the advancing front method (with Laplacian smoothing). Here some of the links I visited: 1.) public domain (TMG seems like a good starting point eventhough it is in C) 2.) general advancing front 3.) thesis 4.) modulef

[URL=http://imageshack.us/photo/my-images/842/fluxlines.jpg/][/URL]

Uploaded with [URL=http://imageshack.us]ImageShack.us[/URL]book

Hello Jjgermis,

Your second link outlines pretty well the method behind the advancing front method. Rather than try and implement someone else's code just go ahead and write your own! It is much more fun and rewarding.

One point in this link, they talk of choosing a target point to create equilateral triangles. That's fine if all your edge lengths are at the desired mesh size. But they are correct in suggesting you start at the smallest length, then in all probability this edge length is smaller than your target mesh size. In which case you want to expand the edges by some pre-determined expansion factor, thus you create an isosceles triangle with the two new sides longer than the base length in an attempt to gradually grow the mesh size.

Laplacian smoothing is a fancy name for one of the most basic and simple routines you can imagine!

Therefore writing this code yourself should not be a problem and of course you have no potential issues with hijacking someone else's code.

Your plot was interesting. I know nothing about electrical flux, but I was intrigued that your contours were either blue or red, why not a rainbow of colours? With rainbow coloured contours you could then use colour filling which would improve visibility.

Hi John

The (direct) calculated quantity is the magnetic vector potential with the unit V.s/m. This ausually has little meaning in this form. The flux lines is actually only to see if the solution is ok.

Therefore I distinguish only between the north pole (red) and south pole (blue). In the air gap this is ideally sinusoidal. With plots like the flux density we use a proper colour scheme.

I am looking forward into programing our own mesh!

PS: distmesh seems like a good starting point for looking into the method. In the dissertation there are good examples as you explained them, i.e. mesh and then apply some Laplacian smoothing.

You guys use a lot of triangles. This sort of field distribution works a treat with isoparametric elements. I'd bet you could do as good a job with a couple of dozen 8-node elements ....

Eddie

Hi Eddie,

you are correct with your comment. I have little experience on the real reasons for using triangles but can comment from our software:

The development started somewhere in the 80's when memory was limited. At that stage a trade off was made between memory and accuracy. To obtain better results with less elements higher oder elemenst were used. At present memory is no problem and coding for triangles is much easier.

Eddie,

I still use a 4-node thin shell element derived from SAP-IV, which is based on a 4 triangle formulation. It may not be the most accurate for well behaved models but it is a very robust element for varying conditions or non-linear applications. Unfortunately, my 8-node isoparametric element does not posses the robust gene and is more difficult for mesh generation.

John

I think that the triangles also lead to easy to handle matrices. Increasing the oder implies automatically an increase in the density. I can imagine that higer order will increase the profile after renumbering. For speed the sparsity pattern should be as 'diagonal' as possible.

The application is also important. For electrical quantities the triangle works just fine. Our mechanical guys requires complete other meshes. Anyway, finding solutions to partial differential equations is a very interessting subject which leads to very challenging solutions 😄 Furthermore, during the 80's this topic seemed to boomed a lot. Most of the literature I have are from that decade!

[URL=http://imageshack.us/photo/my-images/148/gpssy.jpg/][/URL]

Uploaded with [URL=http://imageshack.us]ImageShack.us[/URL]

Hi Eddie,

Meshes composed of triangles are relatively easy to create. As John Campbell said, generating a quad mesh for an arbitrary multi-sided domain is no easy task! I wrote my own paver mesh routine to do this.

Once Jjgermis has written his advancing front routine to generate triangular elements, it should be possible to convert this into a quad mesher. But until you actually try and code this yourself, it is not possible to fully comprehend the complexity of the task !

Also to mesh solid CAD models, despite what some vendors may like you to believe, the only practical approach is to use a tetrahedral mesher, yet more triangles!

It was partly a 'tongue in cheek' comment, but it does seem to have provoked discussion - which is what I intended. There are, indeed, some very good reasons to prefer lots of triangles to rather fewer quadrilaterals, especially where the field variable(s) change rapidly - a particular problem with reduced order gaussian integration. I too have used the 'quad made from 4 triangles' element (although in plane, not shell, form), and if one ever gets into contouring field values in a curvilinear quadrilateral element, it is easier if you divide it into triangles first and then contour in each.

It is some time since I did much with FE, and at the time, I concluded that a spreadsheet was what I wanted to generate those huge lists of coordinates and connections rather than an automatic mesher. However, when I did play with generating bands of elements, I always found the 8-noded quadratic serendipity element awkward, but the equivalent 9-noded Lagrangian element (with a central node) made the array of nodes more regular. I think that the trick is to start with (say) huge quadrilaterals, which are then infilled with one of a number of standard patterns of smaller elements, rather than let the meshing grow semi-organically as in the algorithms being discussed here.

Eddie

,

Hi John

As a first step I worked through the thesis you advised in a previous thread (so I tried a mesh by hand - one needs to do a lot tracking). The author has a compact and easy to follow style. It seems like the following two options could work:

1.) Advancing front (AFM); and 2.) Delaunay.

The code (posted) above calculates the geometric centre and then build the triangles (pencil part). I thought that 'meshing' each triangle could help, but I think this wont work!

What is the motivation for suggesting AFM rather than Delaunay?

[URL=http://imageshack.us/photo/my-images/14/afmr.jpg/][/URL]

Uploaded with [URL=http://imageshack.us]ImageShack.us[/URL]

Jjgermis, this is what an advancing front mesh looks like in a simple square domain.

'What is the motivation for suggesting AFM rather than Delaunay?' - I thought that AFM was simpler to understand, but by all means you can go down the Delauney route. But if you wish to eventually write a quad mesher then you will have to go back to AFM.

Hi John

Thank you for your advice! At present my main interest is producing results in desinging permanent magnet machines. I was rather hoping for a 'quick' solution to get a better mesh or at least one with less effort. Meshing seems to be a very interessing and challenging field which I would like to explore.

For the moment I found the following meshing code (in Fortran). Both have documentation which makes life a little bit easier. The current mesher (mesh above) is very simple and I think that a first step is to improve on that. Here are the links as a result of my search: 1.) delaundo 2.) GEOMPACK2

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. [URL=http://imageshack.us/photo/my-images/62/notch3.jpg/][/URL]

[URL=http://imageshack.us/photo/my-images/853/ex1z.jpg/][/URL]

Uploaded with [URL=http://imageshack.us]ImageShack.us[/URL]

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.

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

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].

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.

[URL=http://imageshack.us/photo/my-images/534/cmos.jpg/][/URL]

[URL=http://imageshack.us/photo/my-images/825/cmosd.jpg/][/URL]

[URL=http://imageshack.us/photo/my-images/703/regione.jpg/][/URL]

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