Simple 2D mesher (in Fortran)

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 😃

[URL=http://imageshack.us/photo/my-images/13/lsup.jpg/][/URL]

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.

  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

[URL=http://imageshack.us/photo/my-images/813/lsup.jpg/][/URL]

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

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:

    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!

[URL=http://imageshack.us/photo/my-images/854/ipmabc.jpg/][/URL]

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

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.

[URL=http://imageshack.us/photo/my-images/821/newmesh.jpg/][/URL]

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

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

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 😄

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!

    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

[URL=http://imageshack.us/photo/my-images/832/ceslot.jpg/][/URL]

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: