bamg mesh

nika

When I have a triangular mesh generated by bamg I have my data stored in a 1 column vector format. So each data point corresponds to vertices. I'm not sure about the order the data stored..Is it every 3 following points make a triangle, starting with an outer triangle?

Thank you

Veronika

mathieumorlighem

Hi Veronika,

your vector should be of size . The column i will have the value of your data at the location of the ith vertex (md.mesh.x(i), md.mesh.y(i)). The triangles are indexed with md.mesh.elements. The ith triangle is made of the following 3 vertices md.mesh.elements(i,:)

Best

Mathieu

nika

I see, but how these triangles are sorted? For instance the first row in md.mesh.elements(i,:) gives me 3 vertices, so the first triangle, but we is it located? starting at the left corner? or centre?

Also, if I want to numerically estimate a derivative in each point (for instance as advection (F_i+1 -F_i)/dx + (F_j+1 -F_j)/dy) it's quite tricky to do it when my data is in this mesh format. Do you think it's easier to interpolate them into a grid then? ( though in this case I'm interpolating into grid and back into mesh, which gives some additional error)

mathieumorlighem

The elements of the meshed are numbered randomly. You can get the coordinates of the center by doing mean(md.mesh.x(md.mesh.elements(i,:))) and same for y. You can also use the following plotmodel options:

plotmodel(md,'data','elementnumering');

plotmodel(md,'data','hilightelements','highlight',); %To highlight elements 1 2 and 3

For the estimate a derivative, the best is to look at what we do in slope.m (we use P1 nodal functions). If you want to estimate the bed slopes with a 2d model, you can do this:


bed = md.geometry.base;
=GetNodalFunctionsCoeff(md.mesh.elements,md.mesh.x,md.mesh.y);
summation=;
bx=(bed(md.mesh.elements).*alpha)*summation;
by=(bed(md.mesh.elements).*beta)*summation;
b=sqrt(bx.^2+by.^2);

nika

Thanks, Mathieu

Is there any reference for this method? I'm not sure I fully understand why the slope is calculated in this way

mathieumorlighem

If you google Triangle shape function, or linear interpolation in 2d, you will find lots of information. The idea is that when we do a linear interpolation on a triangle mesh, your field is decomposed as follows:

bed(x,y) = sum_{i=0}^{numberofvertices} bed_i phi_i(x,y)

where bed_i are the values of the bed on the mesh vertices and phi_i are the nodal or shape functions. Then to get the derivatives, we do

dbed/dx = sum_{i=0}^{numberofvertices} bed_i dphi_i/dx(x,y)

GetNodalFunctionsCoeff will give you the x and y derivatives of phi_i over each element. And then we use a simple averaging. This is really standard.

Mathieu

nika

When saying dbed/dx = sum_{i=0}^{numberofvertices} bed_i dphi_i/dx(x,y)

do you mean dbed/dx(x,y)= sum_{i=0}^{numberofvertices} bed_i dphi_i/dx(x,y)?

mathieumorlighem

yes

nika

Ok, makes sense.Thanks Mathieu, just one more question: As I understand doing it this way the patial derivatives in each point are calculated based on the average ofthe data in each triangle. So if I want to substract derivatives of a paramter from the parameter I should take not its value, but the average in the trangle this parameter belongs to? Like: F_new=F_old(md_last_step.mesh.elements)*summation- dF/dx?

mathieumorlighem

Well, this is up to you... if your mesh resolution is fine enough, you should get similar results

nika

The reason I'm asking is that it confuses me that calculating derivatives I get a matrix that has a size much larger than the number of vertices I have, so I can't multiply for instance two matrices bed_i and dphi_i/dx(x,y) cause they size is different (3254 rows for the first one and 6289 rows for the second one).

mathieumorlighem

Hi Nika,

because this is a linear interpolation, the derivatives are constant by element. So dphi_i/dx is of size md.mesh.numberofelements. If you want to use an L2 projection to change its size to md.mesh.numberofvertices you can use:

dbeddx = averaging(md,dbeddx,0);

Mathieu

nika

Great! works like magic, thanks Mathieu!