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

../trunk-jpl/src/m/geometry/GetAreas.py

@@ +1 @@
+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
+

../trunk-jpl/src/m/mesh/GetNodalFunctionsCoeff.py

@@ +1 @@
+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
+