source: issm/oecreview/Archive/13393-13976/ISSM-13471-13472.diff@ 13980

Index: ../trunk-jpl/src/m/geometry/GetAreas.py
===================================================================
--- ../trunk-jpl/src/m/geometry/GetAreas.py     (revision 0)
+++ ../trunk-jpl/src/m/geometry/GetAreas.py     (revision 13472)
@@ -0,0 +1,52 @@
+import numpy
+
+def GetAreas(index,x,y,z=numpy.array([])):
+       """
+       GETAREAS - compute areas or volumes of elements
+
+          compute areas of triangular elements or volumes 
+          of pentahedrons
+
+          Usage:
+             areas  =GetAreas(index,x,y);
+             volumes=GetAreas(index,x,y,z);
+
+          Examples:
+             areas  =GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y);
+             volumes=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y,md.z);
+       """
+
+       #get number of elements and number of nodes
+       nels=numpy.size(index,axis=0)
+       nods=numpy.size(x)
+
+       #some checks
+       if numpy.size(y)!=nods or (z and numpy.size(z)!=nods):
+               raise TypeError("GetAreas error message: x,y and z do not have the same length.")
+       if numpy.max(index)>nods:
+               raise TypeError("GetAreas error message: index should not have values above %d." % nods)
+       if (not z and numpy.size(index,axis=1)!=3):
+               raise TypeError("GetAreas error message: index should have 3 columns for 2d meshes.")
+       if (z and numpy.size(index,axis=1)!=6):
+               raise TypeError("GetAreas error message: index should have 6 columns for 3d meshes.")
+
+       #initialization
+       areas=numpy.zeros(nels)
+       x1=x[index[:,0].astype(int)-1]
+       x2=x[index[:,1].astype(int)-1]
+       x3=x[index[:,2].astype(int)-1]
+       y1=y[index[:,0].astype(int)-1]
+       y2=y[index[:,1].astype(int)-1]
+       y3=y[index[:,2].astype(int)-1]
+
+       #compute the volume of each element
+       if not z:
+               #compute the surface of the triangle
+               areas=(0.5*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)))
+       else:
+               #V=area(triangle)*1/3(z1+z2+z3)
+               thickness=numpy.mean(z[index[:,3:6].astype(int)-1])-numpy.mean(z[index[:,0:3].astype(int)-1])
+               areas=(0.5*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)))*thickness
+
+       return areas
+
Index: ../trunk-jpl/src/m/mesh/GetNodalFunctionsCoeff.py
===================================================================
--- ../trunk-jpl/src/m/mesh/GetNodalFunctionsCoeff.py   (revision 0)
+++ ../trunk-jpl/src/m/mesh/GetNodalFunctionsCoeff.py   (revision 13472)
@@ -0,0 +1,58 @@
+import numpy
+
+def GetNodalFunctionsCoeff(index,x,y):
+       """
+       GETNODELFUNCTIONSCOEFF - compute nodal functions coefficients
+
+          Compute the coefficients alpha beta and optionaly gamma of
+          2d triangular elements. For each element, the nodal function
+          is defined as:
+          N(x,y)=sum(i=1:3) alpha_i * x + beta_i * y + gamma_i
+
+          Usage:
+             [alpha beta]=GetNodalFunctionsCoeff(index,x,y);
+             [alpha beta gamma]=GetNodalFunctionsCoeff(index,x,y);
+
+          Example:
+             [alpha beta gamma]=GetNodalFunctionsCoeff(md.mesh.elements,md.mesh.x,md.mesh.y);
+       """
+
+       #make columns out of x and y
+       x=x.rshape(-1)
+       y=y.rshape(-1)
+
+       #get nels and nods
+       nels=numpy.size(index,axis=0)
+       nods=numpy.size(x)
+
+       #some checks
+       if numpy.size(y)!=nods:
+               raise TypeError("GetNodalFunctionsCoeff error message: x and y do not have the same length.")
+       if numpy.max(index)>nods:
+               raise TypeError("GetNodalFunctionsCoeff error message: index should not have values above %d." % nods)
+       if numpy.size(index,axis=1)!=3:
+               raise TypeError("GetNodalFunctionsCoeff error message: only 2d meshes supported. index should have 3 columns.")
+
+       #initialize output
+       alpha=numpy.zeros((nels,3))
+       beta=numpy.zeros((nels,3))
+
+       #compute nodal functions coefficients N(x,y)=alpha x + beta y +gamma
+       x1=x[index[:,0].astype(int)-1]
+       x2=x[index[:,1].astype(int)-1]
+       x3=x[index[:,2].astype(int)-1]
+       y1=y[index[:,0].astype(int)-1]
+       y2=y[index[:,1].astype(int)-1]
+       y3=y[index[:,2].astype(int)-1]
+       invdet=1./(x1*(y2-y3)-x2*(y1-y3)+x3*(y1-y2))
+
+       #get alpha and beta
+       alpha=numpy.hstack((invdet*(y2-y3),invdet*(y3-y1),invdet*(y1-y2)))
+       beta =numpy.hstack((invdet*(x3-x2),invdet*(x1-x3),invdet*(x2-x1)))
+
+       #get gamma if requested
+       gamma=numpy.zeros((nels,3))
+       gamma=numpy.hstack((invdet*(x2*y3-x3*y2),invdet*(y1*x3-y3*x1),invdet*(x1*y2-x2*y1)))
+
+       return alpha,beta,gamma
+