fix: update of conceptual models

using numpt functions
Loop3D · Feb 25, 2022 · 948beef · 948beef
1 parent 0b60347
commit 948beef
Showing 1 changed file with 65 additions and 56 deletions.
diff --git a/LoopStructural/modelling/intrusions/geom_conceptual_models.py b/LoopStructural/modelling/intrusions/geom_conceptual_models.py
@@ -4,9 +4,12 @@
 
 
 def ellipse_function(
-    lateral_contact_data, minP=None, maxP=None, minS=None, maxS=None
+    lateral_contact_data, 
+    minP=None, 
+    maxP=None, 
+    minS=None, 
+    maxS=None
 ):
-    import math
 
     if minP == None:
         minP = lateral_contact_data["coord1"].min()
@@ -18,23 +21,26 @@ def ellipse_function(
         maxS = lateral_contact_data["coord2"].max()
 
     a = (maxP - minP) / 2
-    b = (maxS - minS) / 10
+    b = (maxS - minS) / 2
 
     po = minP + (maxP - minP) / 2
 
     p_locations = lateral_contact_data.loc[:, "coord1"].copy().to_numpy()
 
     s = np.zeros([len(p_locations), 2])
 
-    s[np.logical_and(p_locations>minP, p_locations<maxP),0] =  b * np.sqrt(1 - np.power((p_locations[i] - po) / a, 2)) 
-    s[np.logical_and(p_locations>minP, p_locations<maxP),1] =  -b * np.sqrt(1 - np.power((p_locations[i] - po) / a, 2)) 
+    s[np.logical_and(p_locations>minP, p_locations<maxP),0] =  b * np.sqrt(1 - np.power((p_locations[np.logical_and(p_locations>minP, p_locations<maxP)] - po) / a, 2)) 
+    s[np.logical_and(p_locations>minP, p_locations<maxP),1] =  -b * np.sqrt(1 - np.power((p_locations[np.logical_and(p_locations>minP, p_locations<maxP)] - po) / a, 2)) 
 
-
     return s
 
 
 def rectangle_function(
-    lateral_contact_data, minP=None, maxP=None, minS=None, maxS=None
+    lateral_contact_data, 
+    minP=None, 
+    maxP=None, 
+    minS=None, 
+    maxS=None
 ):
     import math
 
@@ -79,52 +85,55 @@ def parallelepiped_function(
 
 
 def obliquecone_function(
-    othercontact_data,
-    mean_growth=None,
-    minP=None,
-    maxP=None,
-    minS=None,
-    maxS=None,
-    vertex=None,
-):
-    import math
-
-    ps_locations = othercontact_data.loc[:, ["coord1", "coord2"]].to_numpy()
-
-    a = (maxP - minP) / 2  # semi-major axis
-    b = (maxS - minS) / 2  # semi-minor axis
-    a2 = pow(a, 2)
-    b2 = pow(b, 2)
-
-    po = minP + a  # p coordinate of ellipsis centre
-    so = minS + b  # s coordinate of ellipsis centre
-
-    alpha = vertex[0]  # p coordinate of vertex
-    beta = vertex[1]  # g coordinate of vertex
-    gamma = vertex[2]  # s coordinate of vertex
-
-    growth = np.zeros([len(ps_locations), 2])  # container for results
-
-    for i in range(len(ps_locations)):
-        p = ps_locations[i, 0]
-        s = ps_locations[i, 1]
-
-        A = alpha - po
-        B = beta * (p - alpha)
-        C = gamma - so
-        D = beta * (s - gamma)
-
-        F = pow(A * b, 2) + pow(C * a, 2) - a2 * b2
-        G = 2 * (B * A * b2 + C * D * a2)
-        H = pow(b * B, 2) + pow(a * D, 2)
-
-        constant_g2 = F
-        constant_g = -2 * F * beta - G
-        constant_1 = F * pow(beta, 2) + G * beta + H
-
-        discriminant = pow(constant_g, 2) - 4 * constant_g2 * constant_1
-
-        growth[i, 0] = -(constant_g + math.sqrt(discriminant)) / (2 * constant_g2)
-        growth[i, 1] = -(constant_g - math.sqrt(discriminant)) / (2 * constant_g2)
-
-    return growth
+    othercontact_data, 
+    mean_growth=None, 
+    minP=None, 
+    maxP=None, 
+    minS=None, 
+    maxS=None, 
+    vertex=None): 
+
+    ps_locations = othercontact_data.loc[:,['coord1','coord2']].to_numpy()
+
+    minP=1.5*minP
+    maxP=1.5*maxP 
+    minS=1.5*minS
+    maxS=1.5*maxS
+
+    a = (maxP-minP)/2 #semi-major axis
+    b = (maxS-minS)/2 #semi-minor axis
+    a2 = pow(a,2)
+    b2 = pow(b,2)
+
+    po = minP + a #p coordinate of ellipsis centre
+    so = minS + b #s coordinate of ellipsis centre
+
+    alpha = vertex[0] #p coordinate of vertex
+    beta = vertex[2] #g coordinate of vertex
+    gamma = vertex[1] #l coordinate of vertex
+
+    growth = np.zeros([len(ps_locations),2]) #container for results
+
+    p = ps_locations[:,0]
+    s = ps_locations[:,1]
+
+    A = alpha - po
+    B = beta*(p[:] - alpha)
+    C = gamma - so
+    D = beta*(s[:] - gamma)
+
+    F = pow(A*b,2) + pow(C*a,2) - a2*b2
+    G = 2*(B*A*b2 + C*D*a2)
+    H = pow(b*B,2) + pow(a*D,2)
+
+    constant_g2 = F
+    constant_g = -2*F*beta - G
+    constant_1 = F*pow(beta,2) + G*beta + H
+
+    discriminant = pow(constant_g,2) - 4*constant_g2*constant_1
+    discriminant[discriminant < 0] = 0
+
+    growth[:,0] = -(constant_g + np.sqrt(discriminant))/(2*constant_g2)
+    growth[:,1] = -(constant_g - np.sqrt(discriminant))/(2*constant_g2)
+
+    return growth