Changeset 26712

Timestamp:

12/06/21 17:57:34 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: fixed interpolation method

File:

• issm/trunk-jpl/src/jl/md/utils.jl (modified) (2 diffs)

issm/trunk-jpl/src/jl/md/utils.jl

@@ -53 +53 @@
         return output
 end# }}}
-function InterpFromMeshToMesh2d(index::Array,x::Vector,y::Vector,data::Vector,xout::Vector,yout::Vector) #{{{
+function InterpFromMeshToMesh2d(index_data::Array,x_data::Vector,y_data::Vector,data::Vector,xout::Vector,yout::Vector,default::Float64=NaN) #{{{
+
+        #Allocate output
+        nods_out = length(xout)
+        data_out = default*ones(nods_out)
+
+        #Interpolation type
+        data_length = size(data,1)
+        nods_data   = length(x_data)
+        nels_data   = size(index_data,1)
+        if(data_length==nods_data)
+                interpolation_type=1;
+        elseif (data_length==nels_data)
+                interpolation_type=2
+        else
+                error("length of vector data not supported yet. It should be of length (number of nodes) or (number of elements)!")
+        end
+        xmin = minimum(xout); xmax = maximum(xout)
+        ymin = minimum(yout); ymax = maximum(yout)
+
+        for i in 1:nels_data
+
+                #skip element if no overlap
+                if (minimum(x_data[index_data[i,:]]) > xmax) continue end
+                if (minimum(y_data[index_data[i,:]]) > ymax) continue end
+                if (maximum(x_data[index_data[i,:]]) < xmin) continue end
+                if (maximum(y_data[index_data[i,:]]) < ymin) continue end
+
+                #get area of the current element (Jacobian = 2 * area)*/
+                #area =x2 * y3 - y2*x3 + x1 * y2 - y1 * x2 + x3 * y1 - y3 * x1;
+                area = (x_data[index_data[i,2]]*y_data[index_data[i,3]]-y_data[index_data[i,2]]*x_data[index_data[i,3]] 
+                                  +  x_data[index_data[i,1]]*y_data[index_data[i,2]]-y_data[index_data[i,1]]*x_data[index_data[i,2]] 
+                                  +  x_data[index_data[i,3]]*y_data[index_data[i,1]]-y_data[index_data[i,3]]*x_data[index_data[i,1]])
+
+                for j in 1:nods_out
+                        #Get first area coordinate = det(x-x3  x2-x3 ; y-y3   y2-y3)/area
+                        area_1=((xout[j]-x_data[index_data[i,3]])*(y_data[index_data[i,2]]-y_data[index_data[i,3]])
+                                         -  (yout[j]-y_data[index_data[i,3]])*(x_data[index_data[i,2]]-x_data[index_data[i,3]]))/area
+                        #Get second area coordinate =det(x1-x3  x-x3 ; y1-y3   y-y3)/area
+                        area_2=((x_data[index_data[i,1]]-x_data[index_data[i,3]])*(yout[j]-y_data[index_data[i,3]])
+                                          - (y_data[index_data[i,1]]-y_data[index_data[i,3]])*(xout[j]-x_data[index_data[i,3]]))/area
+                        #Get third area coordinate = 1-area1-area2
+                        area_3=1-area_1-area_2
+
+                        if (area_1>=0 && area_2>=0 && area_3>=0)
+                                if (interpolation_type==1)
+                                        #nodal interpolation
+                                        data_out[j]=area_1*data[index_data[i,1]]+area_2*data[index_data[i,2]]+area_3*data[index_data[i,3]];
+                                else
+                                        #element interpolation
+                                        data_out[j]=data[i];
+                                end
+                        end
+                end
+        end
+        return data_out
+
+        #OLD STUFF!!! not working...
+        
         #prepare input arrays
         nods = Cint(length(x))
@@ -126 +184 @@
         end
 end #}}}
+function meshgrid(x::Vector, y::Vector)
+    X = [i for i in x, j in 1:length(y)]
+    Y = [j for i in 1:length(x), j in y]
+    return X, Y
+end