Discussion Forum

December 22, 2010 3:10pm UTC

Calculation of normal vectors on surface

Hi,

I need to calculate normal vectors on a 3D geometry surface to calculate boundary conditions from that (it's a rather complex matlab code). I took the postinterp help to get a starting point and tried to calculate the normal vectors in the center of each element. I used postinterp and tried to determine nx, ny and nz, but, whatever I do, I only get NaNs. In the postinterp help, a parameter s was used to do that, this also works fine in 3D, but I am not sure how I would transform certain points to that parameter... Here's my code:

% Create sphere surface geometry
clear fem
fem.geom=geomcoerce('face',sphere3(1,'pos',[0,0,0]));
geomplot(fem,'edgelabels','on'), axis equal
% Generate mesh
fem.mesh = meshinit(fem);
fem.shape = 2;
fem.xmesh = meshextend(fem);
% Get element triangulation
a=get(fem.mesh,'el');
tri=a{3}.elem;
p=fem.mesh.p;
% Calculate center point of each triangle
pn=(p(:,tri(1,:))+p(:,tri(2,:))+p(:,tri(3,:)))/3;
% Calculate normal vectors
nx=postinterp(fem,'nx',pn)
% this gives only NaNs, why?

Any ideas? Also the interpolation in the mesh points does not work,

nx=postinterp(fem,'nx',fem.mesh.p)

gives me the same NaNs...

Nils

Reply  |  Reply with Quote  |  Send private message  |  Report Abuse

December 22, 2010 7:49pm UTC in response to Ivar Kjelberg

Re: Calculation of normal vectors on surface

Well, it would be ideally if I could calculate a normal vector in the center of each element pointing outwards. That would be a on the boundary surface. So: All I want to do is calculate a surface vector somewhere on the boundary surface of a 3D geometry that is not given by a parameter but directly by its coordinates.

I know I have the nx, ny, nz but I just can't get the calculation of it (postinterp, poseval) on an arbitary surface point working. The doc doesn't help, it just gives examples based on a parameterization which I will not be able to use for my final gemeotry.

Reply  |  Reply with Quote  |  Send private message  |  Report Abuse

December 23, 2010 8:41am UTC in response to Ivar Kjelberg

Re: Calculation of normal vectors on surface

My geometry can be pretty curvy and pretty long :-)

As you said - (nx,ny,nz) are defined per element. But how can I get this information per element?

Reply  |  Reply with Quote  |  Send private message  |  Report Abuse

February 9, 2012 6:02am UTC in response to Ivar Kjelberg

Re: Calculation of normal vectors on surface

Calculation of normal vectors on a surface can be tricky. This is how I understand it. Please feel free to correct me.

COMSOL uses "topological" or "surface" coordinates rather than actual coordinates to calculate the normal vectors. In 2D the topological coordinate is one-dimensional (since it is an edge) and in 3D it is 2-dimensional (since it is a surface).

In my opinion, the best way to calculate normal vectors on the surface seems to be a combination of the following commands:

1. flgeomud - Get up-down sub-domain numbering of geometry faces
2. flgeomfd - Get coordinates and derivatives for geometry faces
3. flgeomfn - Get normals from face derivatives

!! Important !!

POSTINTERP command works well for obtaining normal vector information on edges, however returns "NaN" for many face-normal vectors in 3D simulations. Hence I do not recommend POSTINTERP.

Now the philosophy behind calculating normal vectors :

Let number of subdomains be nd. The empty space outside the model is numbered 0.

Commands "flgeomfd" and "flgeomfn" can easily provide us with the face-normal information. But when we talk of normal vector on a surface it becomes mandatory to know how it is directed. For example if domains 1 and 2 share a face number 3, and if the normal vector say [1, 0, 0] on this surface is calculated from the the commands "flgeomfd" and "flgeomfn", we still do not have the information as to whether this normal is directed from domain 1 to domain 2 or domain 2 to domain 1. In other words, the normal vector information provided by these 2 commands lacks 'direction' information.

This is where the command "flgeomud" comes handy. flgeomud is a matrix of size (nd+1) x nf where nf = number of faces. Thus the information from a geometry g obtained using flgeomud looks like this:

>updown = flgeomud (g)

1 1 2 ... ( This row is for "up" domain )
0 0 1 ... ( This row is for "down" domain )

This means :

a) the normal of the face 1 is directed from sub-domain 0 (empty space) to the sub-domain 1
b) the normal of the face 2 is directed from sub-domain 0 (empty space) to the sub-domain 1
c) the normal of the face 3 is directed from sub-domain 1 to the sub-domain 2

Thus the normal vector [1,0,0] of the shared face 3 we obtained earlier points out "from" domain 1 "to" domain 2. Thus it is an "outward" pointing normal for face 1 when domain 1 is concerned and "inward" pointing normal when domain 2 is concerned.

To illustrate this I shall give an example here which presumes that there is only only one subdomain. Thus, we have

empty space : 0
subdomain : 1 = nd

A simple example is a parallelepiped (cuboid) which has six faces. This geometry has faces that are "shared" with the empty space only. Let us think of a situation where there is some heat generated inside this box and that we want to calculate the total heat energy that escapes through all the walls. For this we need to find all the face normal vectors which point "outward" from the box "into" the empty space.

==============================

% Create a simple cubic block

>b=block3(10,20,30);

>clear s
s.objs={b};
s.name={'Block'};
s.tags={'Block'};
fem.draw=struct('s',s);
fem.geom=geomcsg(fem);

% Obtain number of faces (boundaries)

>nb = (geominfo(b, 'out',{'no'},'od',0:3))(3)

nb =

6

% Obtain up-down info

>updown = flgeomud(b)

updown =

1 1 1 0 0 0
0 0 0 1 1 1


% Print normal vector information as obtained by flgeomfd and flgeomfn commands

>for i = 1:nb
[coord,facederiv{i},c] = flgeomfd(b,i,[0.5;0.5]); % Obtain face derivative at mid-surface
normal{i} = flgeomfn(facederiv{i}); % Obtain normal from face-dreivative
end

>normal

normal =

[0 0 1] [1 0 0] [0 1 0] [0 0 1] [1 0 0] [0 1 0]

% The above information lacks direction. Combine with up-down information and
% convert all inward pointing normals into outward pointing normal vectors

>for i = 1:nb
if (updown(1,i)==1) % if inward pointing normal
normal{i}=-normal{i}; % make it outward pointing normal
end
nx(i) = normal{i}(1); % Store in separate vectors (optional)
ny(i) = normal{i}(2);
nz(i) = normal{i}(3);
end

>normal

normal =

[0 0 -1] [-1 0 0] [0 -1 0] [0 0 1] [1 0 0] [0 1 0]

=============

Now that all the outward normal vectors of all the faces have been attained, the total heat energy escaping the box through walls can be easily calculated.

Hope this helps.

Reply  |  Reply with Quote  |  Send private message  |  Report Abuse