From f423d0d91a2d8aea1e0df23e2c009e516d01d55b Mon Sep 17 00:00:00 2001 From: Sabet Seraj <48863473+sseraj@users.noreply.github.com> Date: Mon, 20 Mar 2023 10:40:50 -0400 Subject: [PATCH] Black format update and RST fixes (#23) --- pysurf/__init__.py | 2 +- pysurf/baseClasses.py | 18 --- pysurf/tecplot_interface.py | 15 -- pysurf/tsurf_component.py | 111 ++----------- pysurf/tsurf_tools.py | 264 +++++++++++-------------------- src/complexify.py | 1 - src/f2py/pyf_preprocessor.py | 1 - tests/test_nodal_normals.py | 1 - tests/test_tsurf_component.py | 7 - tests/test_tsurf_curve.py | 2 - tests/test_tsurf_intersection.py | 3 - tests/test_tsurf_remesh.py | 2 - tests/test_tsurf_reorder.py | 2 - 13 files changed, 104 insertions(+), 325 deletions(-) diff --git a/pysurf/__init__.py b/pysurf/__init__.py index 8e9321c..79ca6a0 100644 --- a/pysurf/__init__.py +++ b/pysurf/__init__.py @@ -1,3 +1,3 @@ -__version__ = "1.3.1" +__version__ = "1.3.2" from .tsurf_component import TSurfGeometry, TSurfCurve diff --git a/pysurf/baseClasses.py b/pysurf/baseClasses.py index e8d4f08..89b89fc 100644 --- a/pysurf/baseClasses.py +++ b/pysurf/baseClasses.py @@ -223,7 +223,6 @@ def accumulate_reverseADSeeds(self, coorb=None, curveCoorb=None): # MANIPULATOR INTERFACE METHODS def assign_manipulator(self, GMObj): - """ This function assign a geometry manipulation object (such ad DVGeo) to the current Geometry object. @@ -241,7 +240,6 @@ def assign_manipulator(self, GMObj): # Only the root proc will embed the pySurf nodes into the manipulator if self.myID == 0: - print("") print("Assigning ", self.name, " to manipulator object") @@ -262,7 +260,6 @@ def assign_manipulator(self, GMObj): # Now we need to assign every curve to the manipulator as well for curveName in self.curves: - # Generate name for the curve point set ptSetName = self.name + ":curveNodes:" + curveName @@ -279,12 +276,10 @@ def assign_manipulator(self, GMObj): print("") else: - # The other processors will just save the reference to the manipulator self.manipulator = GMObj def manipulator_addPointSet(self, coor, ptSetName): - """ This method adds an extra point set to the manipulator object. @@ -302,7 +297,6 @@ def manipulator_addPointSet(self, coor, ptSetName): self.manipulator.addPointSet(coor, ptSetName) def manipulator_update(self, ptSetName=None): - """ This method will call the update functions from the manipulator object accordingly to update the geometry based on the new set of design variables. @@ -320,7 +314,6 @@ def manipulator_update(self, ptSetName=None): # Only the root proc should update the triangulated surfaces and curves if self.myID == 0: - # Update surface nodes print("Updating triangulated surface nodes from ", self.name) coor = self.manipulator.update(self.ptSetName) @@ -329,7 +322,6 @@ def manipulator_update(self, ptSetName=None): # Now we need to update every curve as well for curveName in self.curves: - # Update curve nodes print("Updating nodes from curve ", curveName) coor = self.manipulator.update(self.curves[curveName].ptSetName) @@ -338,7 +330,6 @@ def manipulator_update(self, ptSetName=None): # Finally, all procs update the embeded nodes, if they have one if ptSetName is not None: - if self.myID == 0: print("Updating embedded nodes from ", self.name, " under ptSet ", ptSetName) @@ -349,7 +340,6 @@ def manipulator_update(self, ptSetName=None): print("Done") def manipulator_forwardAD(self, xDVd, ptSetName=None): - """ This method uses forward AD to propagate derivative seeds from the design variables to the surface mesh coordinates. @@ -372,7 +362,6 @@ def manipulator_forwardAD(self, xDVd, ptSetName=None): # Now we need to update every curve as well for curveName in self.curves: - # Update curve nodes print("Updating nodes from curve ", curveName) coord = self.manipulator.totalSensitivityProd(xDVd, self.curves[curveName].ptSetName) @@ -382,7 +371,6 @@ def manipulator_forwardAD(self, xDVd, ptSetName=None): # Finally, all procs update the embeded nodes, if they have one if ptSetName is not None: - if self.myID == 0: print("Updating embedded nodes from ", self.name, " under ptSet ", ptSetName) @@ -393,7 +381,6 @@ def manipulator_forwardAD(self, xDVd, ptSetName=None): print("Done") def manipulator_reverseAD(self, xDVb, ptSetName=None, comm=None, clean=True): - """ This method uses reverse AD to propagate derivative seeds from the surface mesh coordinates to the design variables. @@ -404,7 +391,6 @@ def manipulator_reverseAD(self, xDVb, ptSetName=None, comm=None, clean=True): """ if self.myID == 0: - # Print log print("") print("Reverse AD call to manipulator of ", self.name, " object") @@ -420,7 +406,6 @@ def manipulator_reverseAD(self, xDVb, ptSetName=None, comm=None, clean=True): # Now we need to update every curve as well for curveName in self.curves: - # Update curve nodes print("Updating nodes from curve ", curveName) xDVb_curr = self.manipulator.totalSensitivity(curveCoorb[curveName], self.curves[curveName].ptSetName) @@ -429,7 +414,6 @@ def manipulator_reverseAD(self, xDVb, ptSetName=None, comm=None, clean=True): # Finally, all procs update the embeded nodes, if they have one if ptSetName is not None: - if self.myID == 0: print("Updating embedded nodes from ", self.name, " under ptSet ", ptSetName) @@ -440,7 +424,6 @@ def manipulator_reverseAD(self, xDVb, ptSetName=None, comm=None, clean=True): print("Done") def manipulator_getDVs(self): - """ This returns the design variables of the current manipulator. """ @@ -543,7 +526,6 @@ def project(self, xyz): def accumulate_dict(majorDict, minorDict): - """ This function will loop over all keys of minorDict. If majorDict shares the same key, then we will accumulate (add) the values into diff --git a/pysurf/tecplot_interface.py b/pysurf/tecplot_interface.py index 916f079..8fad4a7 100644 --- a/pysurf/tecplot_interface.py +++ b/pysurf/tecplot_interface.py @@ -2,7 +2,6 @@ def write_tecplot_scatter(filename, title, variable_names, data_points): - # Open the data file fid = open(filename, "w") @@ -34,7 +33,6 @@ def write_tecplot_scatter(filename, title, variable_names, data_points): def readTecplotFEdata(fileName): - # Written by Ney Secco # This script will read sections from a Tecplot FE file # FOR NOW THIS JUST READS BAR FE FILES FOR THE CURVE FUNCTIONS @@ -53,7 +51,6 @@ def readTecplotFEdata(fileName): # Loop over the lines to gather data for line in lines: - # Split data data = line.split() @@ -64,7 +61,6 @@ def readTecplotFEdata(fileName): # Detect if we reached a new section definition if data[0].lower() == "zone": - # Gather name of the new zone zoneName = line.split('"')[1] @@ -85,7 +81,6 @@ def readTecplotFEdata(fileName): # Check if we have a data line if len(data) == 3: - # Try to convert the elements of this line to numbers. # If this doesn't work, it means we do not have a data line. try: @@ -95,7 +90,6 @@ def readTecplotFEdata(fileName): # Check if we have a connectivity line if len(data) == 2: - # Try to convert the elements of this line to numbers. # If this doesn't work, it means we do not have a connectivity line. try: @@ -112,7 +106,6 @@ def readTecplotFEdata(fileName): def readTecplotCurves(fileName): - """ This function will read a curve definition from a Tecplot FE data file and assign coor and barsConn. @@ -126,7 +119,6 @@ def readTecplotCurves(fileName): # Create curves for every section curves = [] for secID in range(len(sectionName)): - # Gather data # The -1 is to adjust connectivities to Python indexing, which starts at zero. coor = np.array(sectionData[secID]) @@ -150,7 +142,6 @@ def __init__(self, coor, barsConn, curveName): self.name = curveName def export_tecplot(self, fileName="curve"): - writeTecplotFEdata(self.coor, self.barsConn, self.name, fileName) @@ -159,7 +150,6 @@ def export_tecplot(self, fileName="curve"): def readTecplotFEdataSurf(fileName): - # Written by Ney Secco # This script will read sections from a Tecplot FE file with quads written by pysurf. # The Tecplot file should have a single zone @@ -175,7 +165,6 @@ def readTecplotFEdataSurf(fileName): # Loop over the lines to gather data for line in lines: - # Split data data = line.split() @@ -186,7 +175,6 @@ def readTecplotFEdataSurf(fileName): # Check if we have a coordinate line if len(data) == 3: - # Try to convert the elements of this line to numbers. # If this doesn't work, it means we do not have a data line. try: @@ -196,7 +184,6 @@ def readTecplotFEdataSurf(fileName): # Check if we have a quad connectivity line if len(data) == 4: - # Try to convert the elements of this line to numbers. # If this doesn't work, it means we do not have a connectivity line. try: @@ -226,7 +213,6 @@ def readTecplotFEdataSurf(fileName): def writeTecplotFEdata(coor, barsConn, curveName, fileName): - # This script will write sections from a Tecplot FE file # Written by John Hwang. Adapted by Ney Secco. @@ -273,7 +259,6 @@ def writeTecplotFEdata(coor, barsConn, curveName, fileName): def writeTecplotSurfaceFEData(coor, triaConn, quadsConn, surfName, fileName): - """ This method will export the triangulated surface data into a tecplot format. The will write all elements as quad data. The triangle elements will be exported diff --git a/pysurf/tsurf_component.py b/pysurf/tsurf_component.py index 0148abf..fa03b59 100644 --- a/pysurf/tsurf_component.py +++ b/pysurf/tsurf_component.py @@ -10,6 +10,8 @@ class TSurfGeometry(Geometry): """ + A class for handling a triangulated surface geometry. + Parameters ---------- fileName : string @@ -34,7 +36,6 @@ class TSurfGeometry(Geometry): """ def __init__(self, fileName, sectionsList=None, comm=None, name=None, dtype=float): - # Initialize the base classs super().__init__(comm) @@ -77,12 +78,10 @@ def __init__(self, fileName, sectionsList=None, comm=None, name=None, dtype=floa # Only the root proc will work here if self.myID == 0: - # Get file extension fileExt = os.path.splitext(os.path.basename(fileName))[1] if fileExt == ".cgns": - # Read CGNS file self.coor, sectionDict = tst.getCGNSsections(fileName, self.comm) @@ -98,7 +97,6 @@ def __init__(self, fileName, sectionsList=None, comm=None, name=None, dtype=floa tst.initialize_curves(self, sectionDict, sectionsList, dtype=self.dtype) elif fileExt == ".plt": - # Read Tecplot file self.coor, triaConn, quadsConn = tecplot_interface.readTecplotFEdataSurf(fileName) @@ -122,7 +120,6 @@ def __init__(self, fileName, sectionsList=None, comm=None, name=None, dtype=floa tst.initialize_surface(self) else: - # Other procs create dummy data self.coor = np.zeros((1, 3)) @@ -137,7 +134,6 @@ def rename(self, name): # Only the root proc will work here if self.myID == 0: - # Deallocate previous tree self.adtAPI.adtdeallocateadts(self.name) @@ -164,7 +160,6 @@ def remove_curve(self, name): # Only the root proc will work here if self.myID == 0: - self.curves.pop(name) def rename_curve(self, oldName, newName): @@ -174,21 +169,17 @@ def rename_curve(self, oldName, newName): # Only the root proc will work here if self.myID == 0: - self.curves[newName] = self.curves.pop(oldName) self.curves[newName].rename(newName) def update(self, coor=None): - """ This function updates the nodal coordinates used by the surface object. """ # Only the root proc will work here if self.myID == 0: - if coor is not None: - # First check if we have the same number of new coordinates if not self.coor.shape == coor.shape: print("") @@ -204,30 +195,24 @@ def update(self, coor=None): tst.update_surface(self) def translate(self, x, y, z): - # Only the root proc will work here if self.myID == 0: - tst.translate(self, x, y, z) tst.update_surface(self) for curve in self.curves.values(): curve.translate(x, y, z) def scale(self, factor): - # Only the root proc will work here if self.myID == 0: - tst.scale(self, factor) tst.update_surface(self) for curve in self.curves.values(): curve.scale(factor) def rotate(self, angle, axis, point=None): - # Only the root proc will work here if self.myID == 0: - tst.rotate(self, angle, axis, point) tst.update_surface(self) for curve in self.curves.values(): @@ -237,7 +222,6 @@ def rotate(self, angle, axis, point=None): # SURFACE PROJECTION METHODS def project_on_surface(self, xyz): - """ This function will compute projections and surface Normals @@ -310,38 +294,33 @@ def project_on_surface(self, xyz): return xyzProj, normProj, projDict def project_on_surface_d(self, xyz, xyzd, xyzProj, normProj, projDict): - """ This function will compute derivatives of the projection algorithm in forward mode. Parameters ---------- - xyz: float[numPts, 3] + xyz : float[numPts, 3] Coordinates of the points that should be projected. - xyzd: float[numPts, 3] - Derivative seeds for coordinates of the points - that should be projected. + xyzd : float[numPts, 3] + Derivative seeds for coordinates of the points that should be projected. - xyzProj: float[numPts,3] - Coordinates of the projected points (can be obtained - with the original projection function) + xyzProj : float[numPts,3] + Coordinates of the projected points (can be obtained with the original projection function) - normProj: float[numPts,3] - Surface normal at projected points (can be obtained - with the original projection function) + normProj : float[numPts,3] + Surface normal at projected points (can be obtained with the original projection function) - projDict: dictionary - Dictionary containing intermediate values that are - required by the differentiation routines. (can be obtained - with the original projection function) + projDict : dict + Dictionary containing intermediate values that are required by the differentiation routines + (can be obtained with the original projection function) Returns ------- - xyzProjd: float[numPts,3] + xyzProjd : float[numPts,3] Derivative seeds of the coordinates of the projected points - normProjd: float[numPts,3] + normProjd : float[numPts,3] Derivative seeds of the surface normal at projected points @@ -403,7 +382,6 @@ def project_on_surface_d(self, xyz, xyzd, xyzProj, normProj, projDict): return xyzProjd, normProjd def project_on_surface_b(self, xyz, xyzProj, xyzProjb, normProj, normProjb, projDict): - """ This function will compute derivatives of the projection algorithm in reverse mode. @@ -504,7 +482,6 @@ def project_on_surface_b(self, xyz, xyzProj, xyzProjb, normProj, normProjb, proj # CURVE PROJECTION METHODS def project_on_curve(self, xyz, curveCandidates=None): - """ This function will compute projections and surface Normals @@ -573,7 +550,6 @@ def project_on_curve(self, xyz, curveCandidates=None): # dist2, xyzProj, and normProj for curveName in self.curves: if curveName in curveCandidates: - # Run projection code _, _, _, curveMask = self.curves[curveName].project(xyz, dist2, xyzProj, tanProj, elemIDs) @@ -592,7 +568,6 @@ def project_on_curve(self, xyz, curveCandidates=None): return xyzProj, tanProj, curveProjDict def project_on_curve_d(self, xyz, xyzd, xyzProj, tanProj, curveProjDict): - """ This function will compute projections and surface Normals @@ -634,7 +609,6 @@ def project_on_curve_d(self, xyz, xyzd, xyzProj, tanProj, curveProjDict): # dist2, xyzProj, and normProj for curveName in self.curves: if curveName in curveCandidates: - # Identify points that are projected to the current curve to create a mask curveMask = [0] * numPts for ii in range(numPts): @@ -648,7 +622,6 @@ def project_on_curve_d(self, xyz, xyzd, xyzProj, tanProj, curveProjDict): return xyzProjd, tanProjd def project_on_curve_b(self, xyz, xyzProj, xyzProjb, tanProj, tanProjb, curveProjDict): - """ This function will backpropagate projections derivatives back to curve nodes and unprojected points. @@ -670,9 +643,7 @@ def project_on_curve_b(self, xyz, xyzProj, xyzProjb, tanProj, tanProjb, curvePro # Call inverse_evaluate for each component in the list, so that we can update # dist2, xyzProj, and normProj for curveName in self.curves: - if curveName in curveCandidates: - # Identify points that are projected to the current curve to create a mask curveMask = [0] * numPts for ii in range(numPts): @@ -743,7 +714,6 @@ def intersect_d(self, otherGeometry, intCurve, distTol=1e-7): if (self.name == intCurve.extra_data["parentGeoms"][0]) and ( otherGeometry.name == intCurve.extra_data["parentGeoms"][1] ): - # Print log print("") print("Propagating forward AD seeds for intersection curve:") @@ -751,7 +721,6 @@ def intersect_d(self, otherGeometry, intCurve, distTol=1e-7): print("") else: - raise NameError("Trying to run derivative code with incorrect parents") # Then we can call the AD code defined in tsurf_tools.py @@ -776,7 +745,6 @@ def intersect_b(self, otherGeometry, intCurve, distTol=1e-7, accumulateSeeds=Tru if (self.name == intCurve.extra_data["parentGeoms"][0]) and ( otherGeometry.name == intCurve.extra_data["parentGeoms"][1] ): - # Print log print("") print("Propagating reverse AD seeds for intersection curve:") @@ -784,7 +752,6 @@ def intersect_b(self, otherGeometry, intCurve, distTol=1e-7, accumulateSeeds=Tru print("") else: - raise NameError("Trying to run derivative code with incorrect parents") # Copy the intersection curve seeds to avoid any modification @@ -805,7 +772,6 @@ def intersect_b(self, otherGeometry, intCurve, distTol=1e-7, accumulateSeeds=Tru # POINT METHODS def set_points(self, coor): - """ This will replace the nodal coordinates that define the surface. @@ -818,7 +784,6 @@ def set_points(self, coor): self.update(np.array(coor)) def get_points(self): - """ This will give the nodal coordinates that define the surface. """ @@ -946,7 +911,6 @@ def clean_reverseADSeeds(self): self.curves[curveName].coorb[:, :] = 0.0 def set_randomADSeeds(self, mode="both", fixedSeed=True): - """ This will set random normalized seeds to all variables. This can be used for testing purposes. @@ -961,7 +925,6 @@ def set_randomADSeeds(self, mode="both", fixedSeed=True): # Set forward AD seeds if mode == "forward" or mode == "both": - coord = np.random.random_sample(self.coor.shape) coord = coord / np.sqrt(np.sum(coord**2)) self.coord = coord @@ -973,7 +936,6 @@ def set_randomADSeeds(self, mode="both", fixedSeed=True): # Set reverse AD seeds if mode == "reverse" or mode == "both": - coorb = np.random.random_sample(self.coor.shape) coorb = coorb / np.sqrt(np.sum(coorb**2)) self.coorb = coorb @@ -998,7 +960,6 @@ def set_randomADSeeds(self, mode="both", fixedSeed=True): # VISUALIZATION METHODS def export_tecplot(self, fileName): - """ This method will export the triangulated surface data into a tecplot format. This will write all elements as quad data. The triangle elements will be exported @@ -1009,7 +970,6 @@ def export_tecplot(self, fileName): # Only the root proc will work here if self.myID == 0: - # Call the function that export FE data # we add -1 here to convert from Fortran to Python ordering # The Tecplot function shouldn't assume we have Fortran ordered data. @@ -1050,7 +1010,6 @@ class TSurfCurve(Curve): """ def __init__(self, coor, barsConn, name, mergeTol=1e-7, dtype=float): - # Initialize the base classs super().__init__() @@ -1131,7 +1090,6 @@ def update(self, coor): self.coor = np.array(coor) def flip(self): - self.barsConn = self.barsConn[::-1, ::-1] # Flip the extra data associated with element ordering @@ -1156,7 +1114,6 @@ def rename(self, name): # ================================================================= def project(self, xyz, dist2=None, xyzProj=None, tanProj=None, elemIDs=None): - """ This function will take the points given in xyz and project them to the curve. @@ -1220,7 +1177,6 @@ def project(self, xyz, dist2=None, xyzProj=None, tanProj=None, elemIDs=None): return xyzProj, tanProj, elemIDs, curveMask def project_d(self, xyz, xyzd, xyzProj, xyzProjd, tanProj, tanProjd, elemIDs, curveMask): - """ This function will run forward mode AD to propagate derivatives from inputs (xyz and coor) to outputs (xyzProj). @@ -1278,7 +1234,6 @@ def project_d(self, xyz, xyzd, xyzProj, xyzProjd, tanProj, tanProjd, elemIDs, cu ) def project_b(self, xyz, xyzb, xyzProj, xyzProjb, tanProj, tanProjb, elemIDs, curveMask): - """ This function will run forward mode AD to propagate derivatives from inputs (xyz and coor) to outputs (xyzProj). @@ -1342,7 +1297,6 @@ def project_b(self, xyz, xyzb, xyzProj, xyzProjb, tanProj, tanProjb, elemIDs, cu # ================================================================= def remesh(self, nNewNodes=None, method="linear", spacing="linear", initialSpacing=0.1, finalSpacing=0.1): - """ This function will redistribute the nodes along the curve to get better spacing among the nodes. This assumes that the FE data is ordered. @@ -1431,7 +1385,6 @@ def remesh(self, nNewNodes=None, method="linear", spacing="linear", initialSpaci def remesh_d( self, newCurve, nNewNodes=None, method="linear", spacing="linear", initialSpacing=0.1, finalSpacing=0.1 ): - """ The original curve (before the remesh) should be the one that calls this method. """ @@ -1496,7 +1449,6 @@ def remesh_b( initialSpacing=0.1, finalSpacing=0.1, ): - """ The original curve (before the remesh) should be the one that calls this method. """ @@ -1589,14 +1541,12 @@ def split_d(self, childCurve): # Check if the given curve is actually a child if childCurve.name in list(self.extra_data["splitCurves"].keys()): - # The usedPtsMask of the child curve will help us inherit derivative seeds from # the parent curve usedPtsMask = childCurve.extra_data["usedPtsMask"] for parentNodeID in range(len(usedPtsMask)): - # Get child node that corresponds to the parent node childNodeID = usedPtsMask[parentNodeID] @@ -1614,14 +1564,12 @@ def split_b(self, childCurve): # Check if the given curve is actually a child if childCurve.name in list(self.extra_data["splitCurves"].keys()): - # The usedPtsMask of the child curve will help us inherit derivative seeds from # the parent curve usedPtsMask = childCurve.extra_data["usedPtsMask"] for parentNodeID in range(len(usedPtsMask)): - # Get child node that corresponds to the parent node childNodeID = usedPtsMask[parentNodeID] @@ -1684,16 +1632,13 @@ def merge_d(self, curveDict): # Loop over every curve we want to merge for curveName in self.extra_data["parentCurves"]: - # Check if we have access to this curve if curveName in curveDict: - # Get derivative seeds of the parent curve coord = curveDict[curveName]._get_forwardADSeeds() # Loop over the parent curve nodes to propagate seeds for parentNodeID in range(coord.shape[0]): - # Get corresponding ID in the child nodes childNodeID = linkOld2New[indexOffset + parentNodeID] @@ -1737,13 +1682,10 @@ def merge_b(self, curveDict): # Loop over every curve we want to merge for curveName in self.extra_data["parentCurves"]: - # Check if we have access to this curve if curveName in curveDict: - # Loop over the parent curve nodes to propagate seeds for parentNodeID in range(curveDict[curveName].coor.shape[0]): - # Get corresponding ID in the child nodes childNodeID = linkOld2New[indexOffset + parentNodeID] @@ -1765,13 +1707,10 @@ def merge_b(self, curveDict): # Loop over every curve we want to merge for curveName in self.extra_data["parentCurves"]: - # Check if we have access to this curve if curveName in curveDict: - # Loop over the parent curve nodes to propagate seeds for parentNodeID in range(curveDict[curveName].coor.shape[0]): - # Get corresponding ID in the child nodes childNodeID = linkOld2New[indexOffset + parentNodeID] @@ -1787,7 +1726,6 @@ def merge_b(self, curveDict): # ================================================================= def condense_disconnect_curves(self, guessNode=0): - """ This function takes a curve defined by multiple discontinuous curves (in the FE connectivity sense) and tries to merge them in @@ -1811,11 +1749,10 @@ def condense_disconnect_curves(self, guessNode=0): print("There is no need to condense curves.") print("") else: - # Initialize FEcurves = [] # Create one curve object for each disconnect curve - for (curveID, currConn) in enumerate(sortedConn): + for curveID, currConn in enumerate(sortedConn): curveName = "FEcurve_%03d" % curveID FEcurves.append(TSurfCurve(curveName, self.coor, currConn)) @@ -1831,7 +1768,6 @@ def condense_disconnect_curves(self, guessNode=0): # Loop to chain nodes while len(coorList) > 0: - # Get coordinates of the first and last nodes of the chain firstNode = newCoor[0] lastNode = newCoor[-1] @@ -1844,7 +1780,6 @@ def condense_disconnect_curves(self, guessNode=0): # Find if which end of the chain has the closest new node if np.min(dist2first) < np.min(dist2last): - # The first node of the chain has the closest point. # So we will add a node before the first point @@ -1858,7 +1793,6 @@ def condense_disconnect_curves(self, guessNode=0): newCoor.insert(0, nextNode) else: - # The last node of the chain has the closest point. # So we will add a node after the last point @@ -1887,7 +1821,6 @@ def condense_disconnect_curves(self, guessNode=0): self.barsConn = newBarsConn def shift_end_nodes(self, criteria="maxX", startPoint=None, curveObject=None): - """ This method will shift the finite element ordering of a periodic curve so that the first and last elements @@ -1926,14 +1859,12 @@ def shift_end_nodes(self, criteria="maxX", startPoint=None, curveObject=None): # First we check if the FE data is ordered and periodic for elemID in range(1, nElem): - # Get node indices prevNodeID = barsConn[elemID - 1, 1] currNodeID = barsConn[elemID, 0] # Check if the FE data is ordered if prevNodeID != currNodeID: - # Print warning print("WARNING: Could not shift curve because it has unordered FE data.") print(" Call FEsort first.") @@ -1941,7 +1872,6 @@ def shift_end_nodes(self, criteria="maxX", startPoint=None, curveObject=None): # The first and last nodes should be the same to have a periodic curve. if barsConn[0, 0] != barsConn[-1, 1]: - # Print warning print("WARNING: Could not shift curve because it is not periodic.") return @@ -1952,7 +1882,6 @@ def shift_end_nodes(self, criteria="maxX", startPoint=None, curveObject=None): # reference node (startNodeID) to reorder the FE data if criteria in ["maxX", "maxY", "maxZ"]: - if criteria == "maxX": coorID = 0 elif criteria == "maxY": @@ -1964,7 +1893,6 @@ def shift_end_nodes(self, criteria="maxX", startPoint=None, curveObject=None): startNodeID = np.argmax(coor[:, coorID]) elif criteria in ["minX", "minY", "minZ"]: - if criteria == "minX": coorID = 0 elif criteria == "minY": @@ -1976,7 +1904,6 @@ def shift_end_nodes(self, criteria="maxX", startPoint=None, curveObject=None): startNodeID = np.argmin(coor[:, coorID]) elif criteria == "startPoint": - # Compute the squared distances between the selected startPoint # the coordinates of the curve. deltas = coor - startPoint @@ -1987,7 +1914,6 @@ def shift_end_nodes(self, criteria="maxX", startPoint=None, curveObject=None): startNodeID = np.argmin(dist_2) elif criteria == "curve": - # Find which node is the closest one to the reference curve startNodeID = closest_node(curveObject.coor, coor) @@ -2004,7 +1930,6 @@ def shift_end_nodes(self, criteria="maxX", startPoint=None, curveObject=None): ) def export_tecplot(self, outputName="curve"): - """ This function will export the current curve in tecplot FE format @@ -2022,7 +1947,6 @@ def export_tecplot(self, outputName="curve"): # The last one of the previous bar in barsConn is the first node of the next bar. def get_points(self): - """ We will just return the nodes that define the curve. Use self.remesh first if you want a different set of points. @@ -2045,7 +1969,6 @@ def get_points(self): return pts def set_points(self, pts): - """ We will just set the nodes that define the curve. """ @@ -2080,7 +2003,6 @@ def set_forwardADSeeds(self, coord): self.coord[self.barsConn[ii, 0], :] = coord[ii, :] def get_forwardADSeeds(self): - # Get the number of elements numElems = self.barsConn.shape[0] @@ -2184,7 +2106,6 @@ def clean_reverseADSeeds(self): self.coorb[:, :] = 0.0 def set_randomADSeeds(self, mode="both", fixedSeed=True): - """ This will set random normalized seeds to all variables. This can be used for testing purposes. @@ -2199,14 +2120,12 @@ def set_randomADSeeds(self, mode="both", fixedSeed=True): # Set forward AD seeds if mode == "forward" or mode == "both": - coord = np.random.random_sample(self.coor.shape) coord = coord / np.sqrt(np.sum(coord**2)) self.coord = coord # Set reverse AD seeds if mode == "reverse" or mode == "both": - coorb = np.random.random_sample(self.coor.shape) coorb = coorb / np.sqrt(np.sum(coorb**2)) self.coorb = coorb diff --git a/pysurf/tsurf_tools.py b/pysurf/tsurf_tools.py index 73a8935..5e81d8b 100644 --- a/pysurf/tsurf_tools.py +++ b/pysurf/tsurf_tools.py @@ -14,7 +14,6 @@ def getCGNSsections(inputFile, comm=MPI.COMM_WORLD): - """ This function opens a CGNS file, reads its sections, and returns a dictionary with selected sections. This code also initializes ADT @@ -103,7 +102,6 @@ def getCGNSsections(inputFile, comm=MPI.COMM_WORLD): # Now we split the connectivity arrays according to the sections iSurf = 0 for surf in surfNames: - # Get indices to slice connectivity array. iTriaStart = surfTriaPtr[iSurf] iTriaEnd = surfTriaPtr[iSurf + 1] @@ -125,7 +123,6 @@ def getCGNSsections(inputFile, comm=MPI.COMM_WORLD): iCurve = 0 for curve in curveNames: - # Get indices to slice connectivity array. iBarsStart = curveBarsPtr[iCurve] iBarsEnd = curveBarsPtr[iCurve + 1] @@ -150,7 +147,6 @@ def getCGNSsections(inputFile, comm=MPI.COMM_WORLD): def read_tecplot_curves(fileName, mergeTol=1e-7): - """ This function reads curve described in tecplot format. @@ -167,7 +163,6 @@ def read_tecplot_curves(fileName, mergeTol=1e-7): # Now we convert every Tecplot data curve to a TSurf curve for curveID in range(len(tecCurves)): - # Get current tecplot curve currTecCurve = tecCurves[curveID] @@ -187,7 +182,6 @@ def read_tecplot_curves(fileName, mergeTol=1e-7): def initialize_surface(TSurfGeometry): - """ This function receives a TSurf component and initializes its surface based on the initial connectivity and coordinates @@ -203,7 +197,6 @@ def initialize_surface(TSurfGeometry): def update_surface(TSurfGeometry): - """ This function receives a TSurf component and initializes its surface ADT based on the current connectivity and coordinates. @@ -241,7 +234,6 @@ def update_surface(TSurfGeometry): def initialize_curves(TSurfGeometry, sectionDict, selectedSections, dtype=float): - """ This function initializes all curves given in sectionDict that are shown in selectedSections. @@ -268,11 +260,9 @@ def initialize_curves(TSurfGeometry, sectionDict, selectedSections, dtype=float) # Now we will initialize curve objects for sectionName in sectionDict: - if ("barsConn" in list(sectionDict[sectionName].keys())) and ( sectionName in selectedSections ): # Then we have a curve section that should be stored - # Get data from current section barsConn = sectionDict[sectionName]["barsConn"] @@ -289,7 +279,6 @@ def initialize_curves(TSurfGeometry, sectionDict, selectedSections, dtype=float) def extract_curves_from_surface(TSurfGeometry, feature="sharpness"): - """ This function will extract features from the surface of TSurfGeometry and return Curve objects. @@ -332,7 +321,6 @@ def extract_curves_from_surface(TSurfGeometry, feature="sharpness"): # Create function to reallocate sharedBarInfo if necessary def reallocate_sharedBarInfo(sharedBarInfo): - # Get old size oldSize = sharedBarInfo.shape[0] @@ -370,7 +358,6 @@ def reallocate_sharedBarInfo(sharedBarInfo): # Loop over all trias for triaID in range(nTria): - if np.mod(triaID, 200) == 0: print("Checking tria", triaID + 1, "of", nTria) @@ -389,16 +376,13 @@ def reallocate_sharedBarInfo(sharedBarInfo): # Check each bar for bar in bars: - # Check if the bar is not defined in the dictionary yet if bar not in list(edge2Elem.keys()): - # Then add a new entry for this bar, containing # the first element we found edge2Elem[bar] = [triaID, None] else: - # The entry already exists. This means we already found # the first element that shares this edge and we just need # to add the second one. @@ -421,7 +405,6 @@ def reallocate_sharedBarInfo(sharedBarInfo): # Loop over all quads for quadID in range(nQuads): - if np.mod(quadID, 200) == 0: print("Checking quad", quadID + 1, "of", nQuads) @@ -443,17 +426,14 @@ def reallocate_sharedBarInfo(sharedBarInfo): # Check each bar for bar in bars: - # Check if the bar is not defined in the dictionary yet if bar not in list(edge2Elem.keys()): - # Then add a new entry for this bar, containing # the first element we found. # Remember that quads are flagged with negative IDS. edge2Elem[bar] = [-quadID, None] else: - # The entry already exists. This means we already found # the first element that shares this edge and we just need # to add the second one. @@ -477,7 +457,6 @@ def reallocate_sharedBarInfo(sharedBarInfo): # Now we transfer the remaining edges in edge2Elem to sharedBarInfo. These remaining bars # are at domain boundaries, so they are linked to a single element only for bar in list(edge2Elem.keys()): - # Check if we should increase the array size. if nEdge == sharedBarInfo.shape[0]: sharedBarInfo = reallocate_sharedBarInfo(sharedBarInfo) @@ -514,7 +493,6 @@ def reallocate_sharedBarInfo(sharedBarInfo): # Loop over the bars to extract features for barID in range(sharedBarInfo.shape[0]): - # Gather data from sharedBarInfo node1 = sharedBarInfo[barID, 0] node2 = sharedBarInfo[barID, 1] @@ -540,7 +518,6 @@ def reallocate_sharedBarInfo(sharedBarInfo): # Now we initiliaze curve objects for each disconnect curve we detected for currCurveConn in selectedBarsConn: - # Increment curve counter curveID = curveID + 1 @@ -567,7 +544,6 @@ def reallocate_sharedBarInfo(sharedBarInfo): def merge_surface_sections(sectionDict, selectedSections): - """ This function merges the connectivity data of all surface sections specified in selectedSections. @@ -586,15 +562,12 @@ def merge_surface_sections(sectionDict, selectedSections): # Loop over the selected surfaces to gather connectivities for sectionName in selectedSections: - if sectionName not in list(sectionDict.keys()): - # User wants a section that is not defined. Print a warning message: print("ERROR: Surface", sectionName, "is not defined. Check surface names in your CGNS file.") quit() elif "triaConnF" in list(sectionDict[sectionName].keys()): # Then we have a surface section - # Assign connectivities if triaConnF is None: # Start new connectivities if we have none @@ -606,7 +579,6 @@ def merge_surface_sections(sectionDict, selectedSections): quadsConnF = np.vstack([quadsConnF, sectionDict[sectionName]["quadsConnF"]]) else: - # The user provided a name that is not a surface section print(sectionName, "is not a surface section.") @@ -623,7 +595,6 @@ def merge_surface_sections(sectionDict, selectedSections): def create_curve_from_points(coor, curveName, periodic=False, mergeTol=1e-7, dtype=float): - """ This function will generate a curve object from a list of points provided by the user. @@ -678,7 +649,6 @@ def create_curve_from_points(coor, curveName, periodic=False, mergeTol=1e-7, dty def merge_curves(curveDict, mergedCurveName, curvesToMerge=None): - """ This function will merge curves in the curveDict dictonary whose name is included in curveNames. @@ -701,10 +671,8 @@ def merge_curves(curveDict, mergedCurveName, curvesToMerge=None): # Loop over every curve we want to merge for curveName in curvesToMerge: - # Check if we have access to this curve if curveName in curveDict: - # Get coor and connectivities of this curve. # Remember to apply the offset to the connectivities coor = curveDict[curveName].coor @@ -728,7 +696,6 @@ def merge_curves(curveDict, mergedCurveName, curvesToMerge=None): def split_curves(curveDict, curvesToBeSplit=None, optionsDict={}, criteria="sharpness"): - """ This function will loop over all curves in curveDict and split them based on a given criteria. The available criteria are: @@ -758,9 +725,7 @@ def split_curves(curveDict, curvesToBeSplit=None, optionsDict={}, criteria="shar # Loop over every curve to check for splits for curveName in list(curveDict.keys()): - if curveName in curvesToBeSplit: - # First we remove the curve component from the dictionary curve = curveDict.pop(curveName) @@ -775,7 +740,6 @@ def split_curves(curveDict, curvesToBeSplit=None, optionsDict={}, criteria="shar def split_curve_single(curve, curveName, optionsDict={}, criteria="sharpness"): - """ This function receives a single curve object, splits it according to a criteria, and then return a new dictionary with split curves. @@ -784,26 +748,28 @@ def split_curve_single(curve, curveName, optionsDict={}, criteria="sharpness"): Parameters ---------- - curve: Curve object -> Curve that will be split + curve : Curve object + Curve that will be split - curveName: string -> Name of the original curve. This name will be used - to name the split curves. + curveName: str + Name of the original curve. This name will be used to name the split curves. - optionsDict: dictionary -> Dictionary containing special options for each split - criteria. The following fields are needed: - For 'sharpness' criteria: -'angle': angle (in degrees) used to split the curve. - the default value is 60 deg. - For 'curve' criteria: -'splittingCurves': list of curve objects. We will split - the current curve in the intersection points with other curves. - For 'node' criteria: -'splittingNodes': list of integers. They represent the node IDs - where we will split the curve. + optionsDict: dict + Dictionary containing special options for each split criteria. The following fields are needed: + For 'sharpness' criteria: -'angle': angle (in degrees) used to split the curve. The default value is 60 deg. + For 'curve' criteria: -'splittingCurves': list of curve objects. + We will split the current curve in the intersection points with other curves. + For 'node' criteria: -'splittingNodes': list of integers. + They represent the node IDs where we will split the curve. - criteria: string -> Criteria that will be used to split curves. The options - available for now are: ['sharpness', 'curve', 'node'] + criteria : str + Criteria that will be used to split curves. + The options available for now are: ['sharpness', 'curve', 'node'] Returns ------- - splitCurveDict: dictionary[curve objects] -> Dictionary containing split curves. + splitCurveDict : dict + Dictionary containing split curves. Ney Secco 2016-08 """ @@ -831,7 +797,6 @@ def split_curve_single(curve, curveName, optionsDict={}, criteria="sharpness"): # The next step will depend on the criteria if criteria == "sharpness": - # This criteria will split the curve if its sharpness (defined here # as the change in tangential direction) is beyond a given threshold. @@ -847,7 +812,6 @@ def split_curve_single(curve, curveName, optionsDict={}, criteria="sharpness"): # Loop over the remaining bars to find sharp kinks for elemID in range(1, nElem): - # Compute tangent direction of the current element currTan = coor[barsConn[elemID, 1], :] - coor[barsConn[elemID, 0], :] currTan = currTan / np.linalg.norm(currTan) @@ -857,7 +821,6 @@ def split_curve_single(curve, curveName, optionsDict={}, criteria="sharpness"): # Check if the angle is beyond a certain threshold if angle > sharpAngle: - # Store the current element as a break position breakList.append(elemID) @@ -867,7 +830,6 @@ def split_curve_single(curve, curveName, optionsDict={}, criteria="sharpness"): # If the curve is periodic, we need to check for breaks between the first and last elements. if barsConn[0, 0] == barsConn[-1, 1]: - # Compute tangent direction of the current element prevTan = coor[barsConn[-1, 1], :] - coor[barsConn[-1, 0], :] prevTan = prevTan / np.linalg.norm(prevTan) @@ -881,7 +843,6 @@ def split_curve_single(curve, curveName, optionsDict={}, criteria="sharpness"): # Check if the angle is beyond a certain threshold if angle > sharpAngle: - # We should break the curve at the initial point # so it will not be periodic anymore isPeriodic = False @@ -916,7 +877,6 @@ def split_curve_single(curve, curveName, optionsDict={}, criteria="sharpness"): # CHECK IF BREAKS WERE DETECTED # We can stop this function earlier if we detect no break points if len(breakList) == 0: - # In this case, we just create a dictionary with the original curve # Initialize dictionary that will contain the split curves @@ -958,7 +918,6 @@ def split_curve_single(curve, curveName, optionsDict={}, criteria="sharpness"): # break point. The remaining part of the curve will be determined depending # if the original curve is periodic or not. for splitID in range(nInnercurves): - # For now copy the original set of nodes splitCoor = coor.copy() @@ -983,7 +942,6 @@ def split_curve_single(curve, curveName, optionsDict={}, criteria="sharpness"): # the number of split curves. The periodic case will have minus one curve compared # to a non-periodic one. if isPeriodic: # Curve is periodic, so we need to define one curve - # For now copy the original set of nodes splitCoor = coor.copy() @@ -1003,7 +961,6 @@ def split_curve_single(curve, curveName, optionsDict={}, criteria="sharpness"): curve.extra_data["splitCurves"][splitCurveName] = [breakList[-1], breakList[0]] else: # Curve is not periodic, so we need to define two curves - # CURVE 0 : before the first break point # For now copy the original set of nodes (the unused ones will be removed) @@ -1056,7 +1013,6 @@ def split_curve_single(curve, curveName, optionsDict={}, criteria="sharpness"): def compute_intersections(TSurfGeometryList, distTol=1e-7, comm=MPI.COMM_WORLD): - """ This function will just compute pair-wise intersections for all components given in TSurfGeometryList. @@ -1079,7 +1035,6 @@ def compute_intersections(TSurfGeometryList, distTol=1e-7, comm=MPI.COMM_WORLD): # Call intersection function for each pair for ii in range(numGeometries): for jj in range(ii + 1, numGeometries): - # Compute new intersections for the current pair newIntersections = _compute_pair_intersection(TSurfGeometryList[ii], TSurfGeometryList[jj], distTol, comm) @@ -1098,7 +1053,6 @@ def compute_intersections(TSurfGeometryList, distTol=1e-7, comm=MPI.COMM_WORLD): def _compute_pair_intersection(TSurfGeometryA, TSurfGeometryB, distTol): - """ This function finds intersection curves between components A and B and returns a list of curve objects corresponding to these intersections. @@ -1132,7 +1086,6 @@ def _compute_pair_intersection(TSurfGeometryA, TSurfGeometryB, distTol): # Retrieve results from Fortran if we have an intersection if np.max(arraySizes[1:]) > 0: - # Second Fortran call to retrieve data from the CGNS file. intersectionArrays = TSurfGeometryA.intersectionAPI.retrievedata(*arraySizes) @@ -1146,7 +1099,6 @@ def _compute_pair_intersection(TSurfGeometryA, TSurfGeometryB, distTol): parentTria = np.array(intersectionArrays[2]).T - 1 else: - # Assign empty arrays coor = np.zeros((0, 3), dtype=dtype) barsConn = np.zeros((0, 2)) @@ -1159,7 +1111,6 @@ def _compute_pair_intersection(TSurfGeometryA, TSurfGeometryB, distTol): Intersections = [] if len(coor) > 0: - # Sort FE data. After this step, barsConn may become a list if we # have two disconnect intersection curves. barsConn, newMap = FEsort(barsConn.tolist()) @@ -1174,8 +1125,7 @@ def _compute_pair_intersection(TSurfGeometryA, TSurfGeometryB, distTol): # Set intersection counter intCounter = -1 - for (currConn, currMap) in zip(barsConn, newMap): - + for currConn, currMap in zip(barsConn, newMap): # Increment counter intCounter = intCounter + 1 @@ -1207,7 +1157,6 @@ def _compute_pair_intersection(TSurfGeometryA, TSurfGeometryB, distTol): def _compute_pair_intersection_d(TSurfGeometryA, TSurfGeometryB, intCurve, coorAd, coorBd, distTol): - """ This function is the backward mode of the intersection computation. This can only be called when we already have the intersection curve. @@ -1262,7 +1211,6 @@ def _compute_pair_intersection_d(TSurfGeometryA, TSurfGeometryB, intCurve, coorA def _compute_pair_intersection_b(TSurfGeometryA, TSurfGeometryB, intCurve, coorIntb, distTol): - """ This function is the backward mode of the intersection computation. This can only be called when we already have the intersection curve. @@ -1320,7 +1268,6 @@ def _compute_pair_intersection_b(TSurfGeometryA, TSurfGeometryB, intCurve, coorI def formatStringArray(fortranArray): - """ This functions receives a numpy array of characters coming from Fortran and returns a list of formatted data. @@ -1337,7 +1284,6 @@ def formatStringArray(fortranArray): # Now we format and concatenate the elements for index in range(shape[0]): - # format string newString = "".join(str(fortranArray[:, index], "utf-8")).strip().lower() @@ -1352,7 +1298,6 @@ def formatStringArray(fortranArray): def FEsort(barsConn): - """ This function can be used to sort connectivities coming from the CGNS file @@ -1379,7 +1324,6 @@ def FEsort(barsConn): newMap = [[i] for i in range(nBars)] while keep_searching: - # If nothing happens, we will get out of the loop keep_searching = False @@ -1397,7 +1341,6 @@ def FEsort(barsConn): # NOTE: An "Element" here may be a concatenation of several bar elements while len(oldConn) > 0: - # Pop another element from oldConn oldElement = oldConn.pop(0) oldElemMap = oldMap.pop(0) @@ -1407,10 +1350,8 @@ def FEsort(barsConn): # Now we check if we can fit this element into any new connectivity newElemCounter = 0 - for (newElement, newElemMap) in zip(newConn, newMap): - + for newElement, newElemMap in zip(newConn, newMap): if oldElement[0] == newElement[0]: - # We will flip the old element and place it at the beginning of the # new connectivity newConn[newElemCounter] = oldElement[::-1] + newElement[1:] @@ -1428,7 +1369,6 @@ def FEsort(barsConn): break elif oldElement[0] == newElement[-1]: - # Just append the old element newConn[newElemCounter] = newElement[:-1] + oldElement @@ -1445,7 +1385,6 @@ def FEsort(barsConn): break elif oldElement[-1] == newElement[0]: - # Place the old element at the beginning newConn[newElemCounter] = oldElement + newElement[1:] @@ -1462,7 +1401,6 @@ def FEsort(barsConn): break elif oldElement[-1] == newElement[-1]: - # Flip and append the old element newConn[newElemCounter] = newElement[:-1] + oldElement[::-1] @@ -1482,7 +1420,6 @@ def FEsort(barsConn): newElemCounter = newElemCounter + 1 if not linked_element: # We have an element that is not connected to anyone - # Define a new curve in newConn newConn.append(oldElement) @@ -1499,7 +1436,6 @@ def FEsort(barsConn): # Fill FE data for each curve for curveID in range(len(newConn)): - # Get current curve curve = newConn[curveID] @@ -1508,7 +1444,6 @@ def FEsort(barsConn): # Assign FE connectivity for pointID in range(1, len(curve)): - # Get point indices prevPoint = curve[pointID - 1] currPoint = curve[pointID] @@ -1523,7 +1458,6 @@ def FEsort(barsConn): # Now we do a final check to remove degenerate bars (bars that begin and end at the same point) # Take every disconnect curve in newConnFE: for curveID in range(len(newConnFE)): - # We convert the connectivity array to list so we can 'pop' elements curveFE = newConnFE[curveID].tolist() @@ -1534,13 +1468,11 @@ def FEsort(barsConn): FEcounter = 0 while FEcounter < len(curveFE): - # Get current finite element FE = curveFE[FEcounter] # Check if the start and end points are the same if FE[0] == FE[1]: - # Remove FE curveFE.pop(FEcounter) @@ -1548,7 +1480,6 @@ def FEsort(barsConn): currMap.pop(FEcounter) else: - # Increment counter FEcounter = FEcounter + 1 @@ -1558,16 +1489,13 @@ def FEsort(barsConn): # Now remove empty curves curveID = 0 while curveID < len(newConnFE): - # Check if current curve has no element if len(newConnFE[curveID].tolist()) == 0: - # If this is the case, remove the curve and its mapping newConnFE.pop(curveID) newMap.pop(curveID) else: - # Increment counter curveID = curveID + 1 @@ -1579,7 +1507,6 @@ def FEsort(barsConn): def remove_unused_points(coor, triaConnF=None, quadsConnF=None, barsConn=None): - """ This function will removed unused points from coor and also update all connectivities. @@ -1645,10 +1572,8 @@ def remove_unused_points(coor, triaConnF=None, quadsConnF=None, barsConn=None): # Now we fill the points of the cropped array for pointID in range(nPoints): - # Check if the point is used if usedPtsMask[pointID] == 1: - # Increment counter cropPointID = cropPointID + 1 @@ -1661,7 +1586,6 @@ def remove_unused_points(coor, triaConnF=None, quadsConnF=None, barsConn=None): usedPtsMask[pointID] = cropPointID else: - # Unused points receive a negative flag usedPtsMask[pointID] = -1 @@ -1697,7 +1621,6 @@ def remove_unused_points(coor, triaConnF=None, quadsConnF=None, barsConn=None): def detect_feature(node1, node2, element1, element2, coor, triaConnF, quadsConnF, feature): - """ This function checks if the bar that connects node1 and node2, and is shared between element1 and element2 has a desired feature. This is @@ -1738,7 +1661,6 @@ def detect_feature(node1, node2, element1, element2, coor, triaConnF, quadsConnF # Check the criteria if feature == "sharpness": - # We need to compute the angle between the normals that share the bar # Let's get two inplane vectors for each element first: @@ -1757,7 +1679,6 @@ def detect_feature(node1, node2, element1, element2, coor, triaConnF, quadsConnF return featureIsDetected elif element2 > 0: # We have a tria - # Get element nodes node1Coor = coor[triaConnF[element2, 0] - 1, :] node2Coor = coor[triaConnF[element2, 1] - 1, :] @@ -1768,7 +1689,6 @@ def detect_feature(node1, node2, element1, element2, coor, triaConnF, quadsConnF v22 = node3Coor - node1Coor elif element2 < 0: # We have a quad - # Get element nodes node1Coor = coor[quadsConnF[-element2, 0] - 1, :] node2Coor = coor[quadsConnF[-element2, 1] - 1, :] @@ -1781,7 +1701,6 @@ def detect_feature(node1, node2, element1, element2, coor, triaConnF, quadsConnF # ELEMENT 1 if element1 > 0: # We have a tria - # Get element nodes node1Coor = coor[triaConnF[element1, 0] - 1, :] node2Coor = coor[triaConnF[element1, 1] - 1, :] @@ -1792,7 +1711,6 @@ def detect_feature(node1, node2, element1, element2, coor, triaConnF, quadsConnF v21 = node3Coor - node1Coor elif element1 < 0: # We have a quad - # Get element nodes node1Coor = coor[quadsConnF[-element1, 0] - 1, :] node2Coor = coor[quadsConnF[-element1, 1] - 1, :] @@ -1823,7 +1741,6 @@ def detect_feature(node1, node2, element1, element2, coor, triaConnF, quadsConnF return featureIsDetected elif feature == "open_ends": - # In this case we check if an edge is shared by a single element only, # as this characterizes an open end. @@ -1839,7 +1756,6 @@ def detect_feature(node1, node2, element1, element2, coor, triaConnF, quadsConnF return featureIsDetected else: - print("ERROR: Feature", feature, "cannot be detected as it is not an option.") quit() @@ -1848,9 +1764,9 @@ def detect_feature(node1, node2, element1, element2, coor, triaConnF, quadsConnF # CURVE MESHING FUNCTIONS # =================================== + # Define hyperbolic tangent generator def hypTanDist(Sp1, Sp2, N): - # This is the hyperbolic tangential spacing developed by ICEMCFD # I got the equations from the ICEM manual # This bunching law is coarser at the middle and finer at the ends @@ -1907,9 +1823,9 @@ def findRootb(b): #################################################### #################################################### + # Define tangent generator def tanDist(Sp1, Sp2, N): - # This is the tangential spacing developed by Ney Secco # This bunching law is coarser at the ends and finer at the middle # of the interval, just like shown below: @@ -1979,9 +1895,9 @@ def func(P): #################################################### #################################################### + # Define cubic generator def cubicDist(Sp1, Sp2, N): - # Sp1: initial spacing (within the [0,1] interval) # Sp2: final spacing (within the [0,1] interval) # N: number of nodes @@ -2111,9 +2027,9 @@ def rotate(Geometry, angle, axis, point=None): # NORMALIZATION FUNCTION + # ORIGINAL FUNCTION def normalize(vec): - """ This is a simple function to normalize an array of vectors. Each row of vec will be treated as a vector that should be normalized. @@ -2142,7 +2058,6 @@ def normalize(vec): # FORWARD MODE def normalize_d(vec, vecd): - """ This is the forward differentiation of the normalize routine @@ -2166,7 +2081,6 @@ def normalize_d(vec, vecd): # REVERSE MODE def normalize_b(vec, normalVecb): - """ This is the backward differentiation of the normalize routine @@ -2248,67 +2162,75 @@ def airfoil_intersection( wingGeomName, exportFeatures=False, ): - """ This script works for blunt traling edges - intCurveName: string -> Name of the intersection curve that should be - treated as a wing-body intersection. - This intersection curve should already be in the - manager. manager.intersect() returns curve names. - - LECurve: TSurf curve object -> Curve that follows the wing leading edge. - - UpperTECurve: TSurf curve object -> Curve that follows the wing upper trailing edge. - - LowerTECurve: TSurf curve object -> Curve that follows the wing upper trailing edge. - - mergedCurveName: string -> Name that will be assigned to the intersection curve and - its corresponding collar mesh. - - optionsBodyMesh: Dictionary -> Options to grow the body mesh. The user does not need - to specify boundary conditions. This is done automatically. - Here is an example of dictionary: - optionsBodyMesh = { - 'dStart' : 0.01*2/fact, - 'numLayers' : int(48/2*fact)+1, #Should be /4 - 'extension' : 1.8, - 'epsE0' : 3.5, - 'theta' : -0.5, - 'alphaP0' : 0.25, - 'numSmoothingPasses' : 0, - 'nuArea' : 0.16, - 'numAreaPasses' : 20, - 'sigmaSplay' : 0.3, - 'cMax' : 5.0, - 'ratioGuess' : 10.0, - } - - optionsWingMesh: Dictionary -> Options to grow the wing mesh. The user does not need - to specify boundary conditions. This is done automatically. - Here is an example of dictionary: - - optionsWingMesh = { - 'dStart' : 0.01*2/fact, - 'numLayers' : int(48/2*fact)+1, - 'extension' : 2.2, - 'epsE0' : 12.5, - 'theta' : -0.8, - 'alphaP0' : 0.25, - 'numSmoothingPasses' : 4, - 'nuArea' : 0.16, - 'numAreaPasses' : 0, - 'sigmaSplay' : 0.3, - 'cMax' : 5.0, - 'ratioGuess' : 10.0, - 'guideCurves' : ['curve_te_low','curve_te_upp'], - 'remesh' : True - } - - wingGeomName: string -> Name of the geometry component that represents the wing. - We need this information because pySurf may store either - the wing or the body as the first marching side for the - collar mesh generation. + intCurveName : str + Name of the intersection curve that should be treated as a wing-body intersection. + This intersection curve should already be in the manager. manager.intersect() returns curve names. + + LECurve : TSurf curve object + Curve that follows the wing leading edge. + + UpperTECurve : TSurf curve object + Curve that follows the wing upper trailing edge. + + LowerTECurve : TSurf curve object + Curve that follows the wing upper trailing edge. + + mergedCurveName : string + Name that will be assigned to the intersection curve and its corresponding collar mesh. + + optionsBodyMesh : dict + Options to grow the body mesh. The user does not need to specify boundary conditions. + This is done automatically. Here is an example of dictionary + + .. code-block:: python + + optionsBodyMesh = { + 'dStart' : 0.01*2/fact, + 'numLayers' : int(48/2*fact)+1, #Should be /4 + 'extension' : 1.8, + 'epsE0' : 3.5, + 'theta' : -0.5, + 'alphaP0' : 0.25, + 'numSmoothingPasses' : 0, + 'nuArea' : 0.16, + 'numAreaPasses' : 20, + 'sigmaSplay' : 0.3, + 'cMax' : 5.0, + 'ratioGuess' : 10.0, + } + + + optionsWingMesh : dict + Options to grow the wing mesh. The user does not need to specify boundary conditions. + This is done automatically. Here is an example of dictionary: + + .. code-block:: python + + optionsWingMesh = { + 'dStart' : 0.01*2/fact, + 'numLayers' : int(48/2*fact)+1, + 'extension' : 2.2, + 'epsE0' : 12.5, + 'theta' : -0.8, + 'alphaP0' : 0.25, + 'numSmoothingPasses' : 4, + 'nuArea' : 0.16, + 'numAreaPasses' : 0, + 'sigmaSplay' : 0.3, + 'cMax' : 5.0, + 'ratioGuess' : 10.0, + 'guideCurves' : ['curve_te_low','curve_te_upp'], + 'remesh' : True + } + + wingGeomName : str + Name of the geometry component that represents the wing. + We need this information because pySurf may store either + the wing or the body as the first marching side for the + collar mesh generation. """ @@ -2342,7 +2264,6 @@ def airfoil_intersection( # Loop over every curve after the split for curveName in splitCurveNames: - # Get pointer to the curve object curve = manager.intCurves[curveName] @@ -2360,20 +2281,17 @@ def airfoil_intersection( # Now we can identify and remesh each curve properly for curveName in splitCurveNames: - # Get pointer to the curve object curve = manager.intCurves[curveName] # The trailing edge curve probably has less nodes than the other ones if curve.numNodes == minNumNodes: - # This is the trailing edge curve. # Just apply an uniform spacing optionsDict = {"nNewNodes": int(numTENodes)} TE_curveName = manager.remesh_intCurve(curveName, optionsDict) else: - # We have an upper or lower skin curve. # First let's identify if the curve is defined from # LE to TE or vice-versa @@ -2390,11 +2308,9 @@ def airfoil_intersection( # Now we can determine if we have upper or lower skin if curr_maxCoor < maxCoor: - # We are at the lower skin if LE_to_TE: - optionsDict = { "nNewNodes": numSkinNodes, "spacing": "hypTan", @@ -2405,7 +2321,6 @@ def airfoil_intersection( LS_curveName = manager.remesh_intCurve(curveName, optionsDict) else: - optionsDict = { "nNewNodes": numSkinNodes, "spacing": "hypTan", @@ -2416,11 +2331,9 @@ def airfoil_intersection( LS_curveName = manager.remesh_intCurve(curveName, optionsDict) else: - # We are at the upper skin if LE_to_TE: - optionsDict = { "nNewNodes": numSkinNodes, "spacing": "hypTan", @@ -2431,7 +2344,6 @@ def airfoil_intersection( US_curveName = manager.remesh_intCurve(curveName, optionsDict) else: - optionsDict = { "nNewNodes": numSkinNodes, "spacing": "hypTan", diff --git a/src/complexify.py b/src/complexify.py index d647ec2..2b3ab8d 100644 --- a/src/complexify.py +++ b/src/complexify.py @@ -324,7 +324,6 @@ def join_lines(i, lines): def fix_line(line, implicit_found): - # skip commented lines if patt_comment.search(line) != None: return (line, implicit_found) diff --git a/src/f2py/pyf_preprocessor.py b/src/f2py/pyf_preprocessor.py index 46aaa5a..cbd622b 100644 --- a/src/f2py/pyf_preprocessor.py +++ b/src/f2py/pyf_preprocessor.py @@ -65,7 +65,6 @@ else: # We have a normal line to write: if cur_mode == "both" or cur_mode == mode: - # We now have to check for real/integer types and process: if "kind=inttype" in orig_lines[i]: orig_lines[i] = orig_lines[i].replace("kind=inttype", "kind=4") diff --git a/tests/test_nodal_normals.py b/tests/test_nodal_normals.py index 9a423d4..6fa4591 100644 --- a/tests/test_nodal_normals.py +++ b/tests/test_nodal_normals.py @@ -7,7 +7,6 @@ class TestNodalNormals(unittest.TestCase): - N_PROCS = 2 def test_nodal_normals(self): diff --git a/tests/test_tsurf_component.py b/tests/test_tsurf_component.py index 2ac7ca8..304acde 100644 --- a/tests/test_tsurf_component.py +++ b/tests/test_tsurf_component.py @@ -17,7 +17,6 @@ class TestTSurfProjection(unittest.TestCase): - N_PROCS = 1 def setUp(self): @@ -89,7 +88,6 @@ def setUpDerivTest(self): return xyz, stepSize_FD, stepSize_CS, xyz_d, coor_d, allCoord, xyzProj_b, normProj_b, tanProj_b def computeSurfaceProjections_AD(self, xyz, xyz_d, coor_d, xyzProj_b, normProj_b): - cube = self.cube # Call projection algorithm @@ -109,7 +107,6 @@ def computeSurfaceProjections_AD(self, xyz, xyz_d, coor_d, xyzProj_b, normProj_b return xyzProj, xyzProj_d, normProj, normProj_d, xyz_b, coor_b def computeSurfaceProjections_pert(self, xyz_pert, coor_pert, CS=False): - if CS: cube = self.cube_CS else: @@ -125,7 +122,6 @@ def computeSurfaceProjections_pert(self, xyz_pert, coor_pert, CS=False): return xyzProj_pert, normProj_pert def computeCurveProjections_AD(self, xyz, xyz_d, allCoor_d, xyzProj_b, tanProj_b): - cube = self.cube # Call projection algorithm @@ -148,7 +144,6 @@ def computeCurveProjections_AD(self, xyz, xyz_d, allCoor_d, xyzProj_b, tanProj_b return xyzProj, xyzProj_d, tanProj, tanProj_d, xyz_b, allCoor_b def computeCurveProjections_pert(self, xyz_pert, allCoor_pert, CS=True): - if CS: cube = self.cube_CS else: @@ -165,7 +160,6 @@ def computeCurveProjections_pert(self, xyz_pert, allCoor_pert, CS=True): return xyzProj, tanProj def test_surface_projection_deriv(self): - cube = self.cube xyz, stepSize_FD, stepSize_CS, xyz_d, coor_d, allCoord, xyzProj_b, normProj_b, tanProj_b = self.setUpDerivTest() @@ -206,7 +200,6 @@ def test_surface_projection_deriv(self): np.testing.assert_allclose(normProj_d, normProj_CS, rtol=3e-6, atol=1e-9) def test_curve_projection_deriv(self): - cube = self.cube xyz, stepSize_FD, stepSize_CS, xyz_d, coor_d, allCoord, xyzProj_b, normProj_b, tanProj_b = self.setUpDerivTest() diff --git a/tests/test_tsurf_curve.py b/tests/test_tsurf_curve.py index 354791a..9cc5f20 100644 --- a/tests/test_tsurf_curve.py +++ b/tests/test_tsurf_curve.py @@ -9,11 +9,9 @@ class TestCurveProjection(unittest.TestCase): - N_PROCS = 1 def test_orig_curve_projection(self): - coor = np.array( [[0.0, 0.0, 0.0], [-1, -1, -1], [-1, -1, -1], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0], [1.0, 1.0, 0.0]] ) diff --git a/tests/test_tsurf_intersection.py b/tests/test_tsurf_intersection.py index bd26e3d..50a9c2c 100644 --- a/tests/test_tsurf_intersection.py +++ b/tests/test_tsurf_intersection.py @@ -6,7 +6,6 @@ class TestCurveIntersection(unittest.TestCase): - N_PROCS = 1 def setUp(self): @@ -27,7 +26,6 @@ def setUp(self): self.comp2_CS = pysurf.TSurfGeometry(cylinderFile, comm=comm, dtype=complex) def test_curve_intersection(self): - comp1 = self.comp1 comp2 = self.comp2 @@ -56,7 +54,6 @@ def test_curve_intersection(self): self.assertEqual(len(intersectionSplit), 27) def test_intersection_deriv(self): - # Define distance tolerance distTol = 1e-7 diff --git a/tests/test_tsurf_remesh.py b/tests/test_tsurf_remesh.py index a3965fb..d6016a0 100644 --- a/tests/test_tsurf_remesh.py +++ b/tests/test_tsurf_remesh.py @@ -4,11 +4,9 @@ class TestRemesh(unittest.TestCase): - N_PROCS = 2 def setUp(self): - # Create simple curve and store within the class coor = np.array([[0.0, 0.2, 0.9], [0.1, 0.3, 0.7], [0.5, 0.6, 0.5], [0.8, 0.7, 0.4], [1.0, 0.9, 0.2]]) diff --git a/tests/test_tsurf_reorder.py b/tests/test_tsurf_reorder.py index 90ea693..5deb416 100644 --- a/tests/test_tsurf_reorder.py +++ b/tests/test_tsurf_reorder.py @@ -4,11 +4,9 @@ class TestReorder(unittest.TestCase): - N_PROCS = 2 def test_reorder(self): - # Define number of nodes in the circle nNodes = 5