Theoretical Electrochemistry

February 10, 2012

Calculating Normal Vectors of Boundaries in COMSOL

Filed under: COMSOL — Murali Venkatraman @ 12:45 pm

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.

Advertisements

Leave a Comment »

No comments yet.

RSS feed for comments on this post. TrackBack URI

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Create a free website or blog at WordPress.com.

%d bloggers like this: