Subject: [FEniCS-users] cylindrical mesh

From: Catherine Micek <[email-masked]>
Date: Tue, 11 Aug 2009 05:20:45 -0500

Hi,

I am working on a 3D finite element code and would like to use a  
cylinder mesh.  Since there is no built-in cylinder mesh, I tried  
generating one from a the built-in UnitCube mesh:

# Create mesh
mesh = UnitCube(2,2,2)

# Translate origin from (0.5, 0.5, 0.5) to (0.0, 0.0)
for x in mesh.coordinates():
         x[0] = x[0] - 0.5
         x[1] = x[1] - 0.5
         x[2] = 4.0*(x[2] - 0.5)

# Create boundary mesh
boundary = BoundaryMesh(mesh)

# Stretch to a circle
for x in boundary.coordinates():
         # Circle parameters
         theta = atan2(x[1],x[0])
         r = 2.0
         x[0] = r*cos(theta)
         x[1] = r*sin(theta)
         x[2] = x[2]

# Move boundary values of mesh
mesh.move(boundary)

This code, however, does not generate any interior nodes.  Also, when  
I look at the .xml file, some identical vertices are numbered twice,  
as though they are being interpreted as different vertices.  So my  
questions are these:

1. Is there a way to correct this code so that the mesh that includes  
interior nodes can be generated?
2. I have Matlab code to generate a cylinder; is there an easy way to  
convert that to a form that can be used by fenics?

Any suggestions would be greatly appreciated -- thanks!

Katy

---
Catherine (Katy) A. Micek
Graduate Assistant, U of MN Mathematics Department
http://www.math.umn.edu/~mice0012

---

From: "Garth N. Wells" <[email-masked]>
Date: Tue, 11 Aug 2009 11:25:42 +0100

Catherine Micek wrote:
> Hi,
> 
> I am working on a 3D finite element code and would like to use 
> a cylinder mesh.  Since there is no built-in cylinder mesh, I tried generating one from a the 
> built-in UnitCube mesh:
> 
> # Create mesh
> mesh = UnitCube(2,2,2)
> 
> This code, however, does not generate any interior nodes.  Also, when I 
> look at the .xml file, some identical vertices are numbered twice, as 
> though they are being interpreted as different vertices.  So my 
> questions are these:  
> 

Not without some work. Also, the quality may be poor when mapping from a 
cube.

> convert that to a form that can be used by fenics?

No. You could modify the Python script dolfin-convert to convert from 
the Matlab format to the DOLFIN XML format.

The simplest approach would be to use a mesh generator that produces a 
format supported by dolfin-convert. Try gmsh

     http://www.geuz.org/gmsh/

Garth

> _http://www.math.umn.edu/~mice0012___
> 
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> FEniCS-users mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/fenics-users

---

From: Johan Hake <[email-masked]>
Date: Tue, 11 Aug 2009 12:43:30 +0200

On Tuesday 11 August 2009 12:25:42 Garth N. Wells wrote:
> > a cylinder mesh.  Since there is no built-in cylinder mesh, I tried

If your matlab code creates all the interior vertices and assign the vertices 
each cell has, you might try to build your mesh using the MeshEditor class in 
PyDOLFIN. I have not used it from python but it "should just work"

  m = Mesh()
  me = MeshEditor()
  me.open(m,3,3)
  me.init_vertices(num_vertices)
  me.add_vertex(0,x0,y0,z0)
  ...
  me.init_cells(num_cells)
  me.add_cell(0,v0,v1,v2,v3)
  ...
  me.close()

  f = File("My_mesh.xml")
  f << m

Johan

---

From: Anders Logg <[email-masked]>
Date: Tue, 11 Aug 2009 12:47:29 +0200

> the Matlab format to the DOLFIN XML format.

There is actually an old MATLAB script in misc/utils/matlab/xmlmesh.m
which writes a mesh to DOLFIN XML format. It hasn't been tested for a
while and probably needs updating to the new DOLFIN XML format. It's
fairly short and simple so you might be able to make the required
changes. Post the updated script if you make it work so we can replace
the old one.

>      http://www.geuz.org/gmsh/

I'd suggest using Netgen. It has something called constructive solid
geometry (CSG) which lets you define meshes by adding and subtracting
simple shapes. Here's a relevant example from the netgen-doc packages:

#
## a cylinder
#
algebraic3d

# cut cylinder by planes:

solid fincyl = cylinder ( 3, 0, 0; -1, 0, 0; 0.5 )
      and plane (0, 0, 0; -1, 0, 0)
      and plane (2, 0, 0; 1, 0, 0);

tlo fincyl;

Unfortunately, Netgen does not support CSG in 2D.

--
Anders

---

From: Marie Rognes <[email-masked]>
Date: Tue, 11 Aug 2009 21:47:34 +0200

> changes. Post the updated script if you make it work so we can replace

I updated this a while back for my own purposes, attached. (Only tested 
in 2d though)

--
Marie

function xmlmesh(filename,p,t)

% XMLMESH - SAVE MATLAB 2D (AND FEMLAB 3D) MESH TO DOLFIN XML FORMAT
%
% Usage: xmlmesh(filename,p,t)
%
%   p - points    (exported from PDE Toolbox)
%   t - triangles (exported from PDE Toolbox)
%
% Copyright (C) 2004 Erik Svensson.
% Licensed under the GNU LGPL Version 2.1.
%
% Modified by Anders Logg 2004-2005.
% Modified by Marie Rognes 2009.
%
% First added:  2004-02-10
% Last changed: 2009-08-11

% Open file
fp = fopen(filename,'w');

np = size(p,2);
nt = size(t,2);

% Write header
fprintf(fp,'<?xml version="1.0" encoding="UTF-8"?>\n\n');
fprintf(fp,'<dolfin xmlns:dolfin="http://www.fenics.org/dolfin">\n');

% 2D mesh
if (size(p,1) == 2)

  % Write nodes
  fprintf(fp,'  <mesh celltype="triangle" dim="2">\n');
  disp('Writing vertices...')
  fprintf(fp,'    <vertices size="%d">\n',np);
  for n=1:np
    fprintf(fp,'      <vertex index="%d" x="%f" y="%f" z="0.0"/>\n', ...
	    n-1, p(1,n), p(2,n));
  end
  fprintf(fp,'    </vertices>\n');

  % Write cells
  disp('Writing cells...')
  fprintf(fp,'    <cells size="%d">\n',nt);
  for n=1:nt
    fprintf(fp,'      <triangle index="%d" v0="%d" v1="%d" v2="%d"/>\n', ...
	    n-1,t(1,n)-1,t(2,n)-1,t(3,n)-1);
  end
  fprintf(fp,'    </cells>\n');
  fprintf(fp,'  </mesh>\n');
  fprintf(fp,'</dolfin>\n');

% 3D mesh
elseif (size(p,1) == 3)

  % Write nodes
  disp('Writing nodes...')
  fprintf(fp,'  <mesh>\n');
  fprintf(fp,'    <vertices size="%d">\n',np);
  for n=1:np
    fprintf(fp,'      <vertex name="%d" x="%f" y="%f" z="%f"/>\n', ...
            n-1,p(1,n),p(2,n),p(3,n));
  end
  fprintf(fp,'    </vertices>\n');

  % Write cells
  disp('Writing cells...')
  fprintf(fp,'    <cells size="%d">\n',nt);
  for n=1:nt
    fprintf(fp,'      <tetrahedron name="%d" n0="%d" n1="%d" n2="%d" n3="%d"/>\n', ...
            n-1,t(1,n)-1,t(2,n)-1,t(3,n)-1,t(4,n)-1);
  end
  fprintf(fp,'    </cells>\n');
  fprintf(fp,'  </mesh>\n');
  fprintf(fp,'</dolfin>\n');

end

% Close file
fclose(fp);
disp('Done')

---

From: Anders Logg <[email-masked]>
Date: Tue, 11 Aug 2009 12:53:26 +0200

On Tue, Aug 11, 2009 at 09:47:34PM +0200, Marie Rognes wrote:
> I updated this a while back for my own purposes, attached. (Only
> tested in 2d though)

Thanks. I've pushed it.

--
Anders

---

From: Catherine Micek <[email-masked]>
Date: Tue, 11 Aug 2009 08:41:50 -0500

>> I updated this a while back for my own purposes, attached. (Only

Ok, thanks for all the suggestions!  I'll look into them.

Katy

>>     fprintf(fp,'      <triangle index="%d" v0="%d" v1="%d" v2="%d"/ 
>> >\n', ...