diff --git a/optvl/om_wrapper.py b/optvl/om_wrapper.py index b75f217..5a4d0dd 100644 --- a/optvl/om_wrapper.py +++ b/optvl/om_wrapper.py @@ -121,7 +121,16 @@ def add_ovl_geom_vars(self, ovl, add_as="inputs", include_airfoil_geom=False): elif add_as == "outputs": self.add_output(geom_key, val=surf_data[surf][key], tags="geom") +def add_ovl_mesh_out_as_output(self, ovl): + surf_data = ovl.get_surface_params() + + meshes,_ = ovl.get_cp_data() + for surf in surf_data: + idx_surf = ovl.surface_names.index(surf) + out_name = f"{surf}:mesh" + self.add_output(out_name, val=meshes[idx_surf], tags="geom_mesh") + def add_ovl_conditions_as_inputs(sys, ovl): # TODO: add all the condition constraints @@ -697,14 +706,18 @@ class OVLMeshReader(om.ExplicitComponent): def initialize(self): self.options.declare("geom_file", types=str) self.options.declare("mass_file", default=None) + self.options.declare("mesh_output",default=False) def setup(self): geom_file = self.options["geom_file"] mass_file = self.options["mass_file"] + mesh_output = self.options["mesh_output"] avl = OVLSolver(geo_file=geom_file, mass_file=mass_file, debug=False) add_ovl_geom_vars(self, avl, add_as="outputs", include_airfoil_geom=True) + if mesh_output: + add_ovl_mesh_out_as_output(self,avl) class Differencer(om.ExplicitComponent): def setup(self): diff --git a/optvl/optvl_class.py b/optvl/optvl_class.py index 357772e..02a765d 100644 --- a/optvl/optvl_class.py +++ b/optvl/optvl_class.py @@ -274,7 +274,7 @@ def __init__( self.avl.loadgeo("") else: raise ValueError("neither a geometry file nor an input options dictionary was specified") - + # todo store the default dict somewhere else # the control surface contraints get added to this array in the __init__ self.conval_idx_dict = { @@ -644,6 +644,9 @@ def check_type(key, avl_vars, given_val): val = optional_header_defaults[key] else: raise ValueError(f"Key {key} not found in input dictionary but is required") + elif key == "title": + # We need to apply this function to the title string so that the tecplot file writing works correctly + val = self._str_to_fort_str(input_dict[key],num_max_char=120) else: val = input_dict[key] @@ -665,11 +668,11 @@ def check_type(key, avl_vars, given_val): # set the gloabl design variable options ndesign = len(input_dict.get("gname", [])) - self.set_avl_fort_arr("CASE_I","NDESIGN", ncontrol) + self.set_avl_fort_arr("CASE_I","NDESIGN", ndesign) if ndesign > self.NGMAX: raise RuntimeError(f"Number of specified design variables exceeds {self.NGMAX}. Raise NGMAX!") - for k in range(ncontrol): + for k in range(ndesign): self.avl.CASE_C.GNAME[k] = input_dict["gname"][k] @@ -680,7 +683,14 @@ def check_type(key, avl_vars, given_val): else: raise RuntimeError(f"Number of specified surfaces, {num_surfs}, exceeds {self.NFMAX}. Raise NFMAX!") - + # Class variable to store the starting index of all meshes. Set to 0 for no mesh. + # We will insert entries into it for duplicate surfaces later but right now it's only for unique surfaces + self.mesh_idx_first = np.zeros(self.get_num_surfaces(),dtype=np.int32) + + # Class variable to store the yoffset for duplicated meshes. This information is needed for correcting retreiving + # a duplicated mesh + self.y_offsets = np.zeros(self.get_num_surfaces(),dtype=np.float64) + # Load surfaces if num_surfs > 0: surf_names = list(input_dict["surfaces"].keys()) @@ -728,15 +738,30 @@ def check_type(key, avl_vars, given_val): "scale": np.array([1.,1.,1.]), "translate": np.array([0.,0.,0.]), "angle": 0.0, + "aincs": np.zeros(num_secs, dtype=np.float64), "wake": True, "albe": True, "load": True, "clcd": np.zeros(6, dtype=np.float64), - "nspans": np.zeros(num_secs, dtype=int), - "sspaces": np.zeros(num_secs, dtype=int), + "nspans": np.zeros(num_secs, dtype=np.int32), + "sspaces": np.zeros(num_secs, dtype=np.float64), "clcdsec": np.zeros((num_secs,6)), "claf": np.ones(num_secs), } + + + ignore_if_mesh = { + "xles", + "yles", + "zles", + "chords", + "nchordwise", + "cspace", + "sspaces" + "sspace", + "nspan", + } + # fmt: on # set some flags based on the options used for this surface @@ -762,11 +787,13 @@ def check_type(key, avl_vars, given_val): ): if key not in surf_dict: - if key in optional_surface_defaults: + if (key == "yduplicate") or (("mesh" in surf_dict) and (key in ignore_if_mesh)): + continue + elif key in optional_surface_defaults: val = optional_surface_defaults[key] else: raise ValueError(f"Key {key} not found in surface dictionary, {surf_name}, but is required") - else: + else: val = surf_dict[key] check_type(key, avl_vars, val) @@ -800,7 +827,7 @@ def check_type(key, avl_vars, given_val): # Load airfoil data sections # Load the Airfoil Section into AVL for j in range(num_secs): - xfminmax = xfminmax_arr[j] + xfminmax = xfminmax_arr[j] # Manually Specify Coordiantes (no camberline verification, only use if you know what you're doing) if "xasec" in surf_dict.keys(): @@ -858,26 +885,27 @@ def check_type(key, avl_vars, given_val): self.set_avl_fort_arr("SURF_GEOM_R", "ZLASEC", np.array([0.0, 0.0]), slicer=slicer_airfoil_flat) self.set_avl_fort_arr("SURF_GEOM_R", "ZUASEC", np.array([0.0, 0.0]), slicer=slicer_airfoil_flat) self.set_avl_fort_arr("SURF_GEOM_R", "CASEC", np.array([0.0, 0.0]), slicer=slicer_airfoil_flat) - - - # --- setup control variables for each section --- # Load control surfaces if "icontd" in surf_dict.keys(): - for j in range(num_secs): + for idx_sec in range(num_secs): # check to make sure this section has control vars - if surf_dict["num_controls"][j] == 0: + if surf_dict["num_controls"][idx_sec] == 0: continue - + for key, avl_vars in self.con_surf_to_fort_var[surf_name].items(): avl_vars_secs = self.con_surf_to_fort_var[surf_name][key] - avl_vars = (avl_vars_secs[0], avl_vars_secs[1], avl_vars_secs[2][j]) + avl_vars = (avl_vars_secs[0], avl_vars_secs[1], avl_vars_secs[2][idx_sec]) if key not in surf_dict: raise ValueError(f"Key {key} not found in surf dictionary, `{surf_name}` but is required") - else: - val = surf_dict[key][j] + else: + # This has to be incremented by 1 for Fortran indexing + if key == "icontd": + val = surf_dict[key][idx_sec] + 1 + else: + val = surf_dict[key][idx_sec] check_type(key, avl_vars, val) self.set_avl_fort_arr(avl_vars[0], avl_vars[1], val, slicer=avl_vars[2]) @@ -887,19 +915,23 @@ def check_type(key, avl_vars, given_val): # --- setup design variables for each section --- # Load design variables if "idestd" in surf_dict.keys(): - for j in range(num_secs): + for idx_sec in range(num_secs): # check to make sure this section has control vars - if surf_dict["num_design_vars"][j] == 0: + if surf_dict["num_design_vars"][idx_sec] == 0: continue for key, avl_vars in self.des_var_to_fort_var[surf_name].items(): avl_vars_secs = self.des_var_to_fort_var[surf_name][key] - avl_vars = (avl_vars_secs[0], avl_vars_secs[1], avl_vars_secs[2][j]) + avl_vars = (avl_vars_secs[0], avl_vars_secs[1], avl_vars_secs[2][idx_sec]) if key not in surf_dict: raise ValueError(f"Key {key} not found in surf dictionary, `{surf_name}` but is required") - else: - val = surf_dict[key][j] + else: + # This has to be incremented by 1 for Fortran indexing + if key == "idestd": + val = surf_dict[key][idx_sec] + 1 + else: + val = surf_dict[key][idx_sec] check_type(key, avl_vars, val) self.set_avl_fort_arr(avl_vars[0], avl_vars[1], val, slicer=avl_vars[2]) @@ -909,10 +941,25 @@ def check_type(key, avl_vars, given_val): # Make the surface if self.debug: print(f"Building surface: {surf_name}") - self.avl.makesurf(idx_surf + 1) # +1 to convert to 1 based indexing + # Load the mesh and make if one is specified otherwise just make + if "mesh" in surf_dict.keys(): + # Check if we have to define the sections for the user + if "iptloc" not in surf_dict.keys(): + surf_dict["iptloc"] = np.zeros(surf_dict["num_sections"],dtype=np.int32) + self.avl.adjust_mesh_spacing(idx_surf+1,surf_dict["mesh"].transpose((2, 0, 1)),surf_dict["iptloc"]) + surf_dict["iptloc"] = surf_dict["iptloc"] - 1 + if "flatten mesh" not in surf_dict.keys(): + surf_dict["flatten mesh"] = True + self.set_mesh(idx_surf, surf_dict["mesh"],surf_dict["iptloc"],flatten=surf_dict["flatten mesh"],update_nvs=True,update_nvc=True) # set_mesh handles the Fortran indexing and ordering + self.avl.makesurf_mesh(idx_surf + 1) #+1 for Fortran indexing + else: + self.avl.makesurf(idx_surf + 1) # +1 to convert to 1 based indexing if "yduplicate" in surf_dict.keys(): self.avl.sdupl(idx_surf + 1, surf_dict["yduplicate"], "YDUP") + # Insert duplicate into the mesh first index array + self.mesh_idx_first = np.insert(self.mesh_idx_first,idx_surf+1,self.mesh_idx_first[idx_surf]) + self.y_offsets[idx_surf] = surf_dict["yduplicate"] self.avl.CASE_I.NSURF += 1 idx_surf += 1 @@ -1023,6 +1070,134 @@ def check_type(key, avl_vars, given_val): # Tell AVL that geometry exists now and is ready for analysis self.avl.CASE_L.LGEO = True + def set_mesh(self, idx_surf: int, mesh: np.ndarray, iptloc: np.ndarray, flatten:bool=True, update_nvs: bool=False, update_nvc: bool=False): + """Sets a mesh directly into OptVL. Requires an iptloc vector to define the indices where the sections are defined. + This is required for many of AVL's features like control surfaces to work properly. OptVL's input routine has multiple + ways of automatically computing this vector. Alternatively, calling the adjust_mesh_spacing subroutine in the Fortran layer + can automatically compute the iptloc vector for a given mesh and number of sections. NOTE: the iptloc input is not differentiated. + Additionally, the length of iptloc cannot change (i.e the number sections cannot change for a surface that's already loaded). + + Args: + idx_surf (int): the surface to apply the mesh to + mesh (np.ndarray): XYZ mesh array (nx,ny,3) + iptloc (np.ndarray): Vector containing the spanwise indicies where each section is defined (num_sections,) + flatten (bool): Should OptVL flatten the mesh when placing vorticies and control points + update_nvs (bool): Should OptVL update the number of spanwise elements for the given mesh + update_nvc (bool): Should OptVL update the number of chordwise elements for the given mesh + """ + nx = copy.deepcopy(mesh.shape[0]) + ny = copy.deepcopy(mesh.shape[1]) + + if len(iptloc) > self.NSMAX: + raise RuntimeError("Length of iptloc cannot exceed NSMAX. Raise NSMAX") + + if update_nvs: + self.avl.SURF_GEOM_I.NVS[idx_surf] = ny-1 + + if update_nvc: + self.avl.SURF_GEOM_I.NVC[idx_surf] = nx-1 + + if flatten: + self.avl.SURF_MESH_L.LMESHFLAT[idx_surf] = True + else: + self.avl.SURF_MESH_L.LMESHFLAT[idx_surf] = False + + # Only add +1 for Fortran indexing if we are not explictly telling the routine to use + # nspans by passing in all zeros + if not (iptloc == 0).all(): + iptloc += 1 + # set iptloc + self.set_avl_fort_arr("SURF_MESH_I","IPTSEC",iptloc,slicer=(idx_surf,slice(None,len(iptloc)))) + + # Compute and set the mesh starting index + if idx_surf != 0: + self.mesh_idx_first[idx_surf] = self.mesh_idx_first[idx_surf-1] + 3*(self.avl.SURF_GEOM_I.NVS[idx_surf-1]+1)*(self.avl.SURF_GEOM_I.NVC[idx_surf-1]+1) + + self.set_avl_fort_arr("SURF_MESH_I","MFRST",self.mesh_idx_first[idx_surf]+1,slicer=idx_surf) + + # Reshape the mesh + # mesh = mesh.ravel(order="C").reshape((3,mesh.shape[0]*mesh.shape[1]), order="F") + mesh = mesh.transpose((1,0,2)).reshape((mesh.shape[0]*mesh.shape[1],3)) + + # Set the mesh + self.set_avl_fort_arr("SURF_MESH_R","MSHBLK",mesh, slicer=(slice(self.mesh_idx_first[idx_surf],self.mesh_idx_first[idx_surf]+nx*ny),slice(0,3))) + + # Flag surface as using mesh geometry + self.avl.SURF_MESH_L.LSURFMSH[idx_surf] = True + + def get_mesh(self, idx_surf: int, get_full_mesh: bool = False, get_iptloc: bool = False): + """ + + Args: + idx_surf (int): the surface to get the mesh for + get_full_mesh (bool): concatenates and returns the meshes for idx_surf and idx_surf + 1, for use with duplicated surfaces + get_iptloc (bool) : should the iptloc vector for the surface be returned + """ + + # Check if surface is using mesh geometry + if not self.avl.SURF_MESH_L.LSURFMSH[idx_surf]: + raise RuntimeError(f"Surface {idx_surf} does not have a mesh assigned!") + + # Get the mesh size + nx = self.avl.SURF_GEOM_I.NVC[idx_surf] + 1 + ny = self.avl.SURF_GEOM_I.NVS[idx_surf] + 1 + + + # Get the mesh + mesh = self.get_avl_fort_arr("SURF_MESH_R","MSHBLK",slicer=(slice(self.mesh_idx_first[idx_surf],self.mesh_idx_first[idx_surf]+nx*ny),slice(0,3))) + mesh = copy.deepcopy(mesh.reshape((ny,nx,3)).transpose((1,0,2))) + # mesh = mesh.transpose((1,0,2)) + + # Check if duplicated + imags = self.get_avl_fort_arr("SURF_I", "IMAGS") + if imags[idx_surf] < 0: + mesh[:,:,1] = -mesh[:,:,1] + self.y_offsets[idx_surf-1] + + # Concatenate with duplicate + if get_full_mesh: + if imags[idx_surf] < 0: + raise RuntimeError(f"Concatenating a duplicated surface, {idx_surf}, with the next surface!") + elif imags[idx_surf+1] > 0: + raise RuntimeError(f"Concatenating with a non duplicated surface, {idx_surf+1}!") + # Get the next mesh + mesh_dup = self.get_avl_fort_arr("SURF_MESH_R","MSHBLK",slicer=(slice(self.mesh_idx_first[idx_surf+1],self.mesh_idx_first[idx_surf+1]+nx*ny),slice(0,3))) + mesh_dup = mesh_dup.reshape((ny,nx,3)).transpose((1,0,2)) + mesh_dup[:,:,1] = -mesh_dup[:,:,1] + self.y_offsets[idx_surf-1] + + # Concatenate them + mesh = np.hstack([mesh,mesh_dup]) + + # Get iptloc + iptloc = None + if get_iptloc: + iptloc = self.get_avl_fort_arr("SURF_MESH_I","IPTSEC",slicer=(idx_surf,slice(None,self.get_num_sections(self.get_surface_names()[idx_surf])))) - 1 + return mesh, iptloc + else: + return mesh + + + # Only add +1 for Fortran indexing if we are not explictly telling the routine to use + # nspans by passing in all zeros + # if not (iptloc == 0).all(): + # iptloc += 1 + # # set iptloc + # self.set_avl_fort_arr("SURF_MESH_I","IPTSEC",iptloc,slicer=(idx_surf,slice(None,len(iptloc)))) + + # Compute and set the mesh starting index + # if idx_surf != 0: + # self.mesh_idx_first[idx_surf] = self.mesh_idx_first[idx_surf-1] + 3*(self.avl.SURF_GEOM_I.NVS[idx_surf-1]+1)*(self.avl.SURF_GEOM_I.NVC[idx_surf-1]+1) + + # self.set_avl_fort_arr("SURF_MESH_I","MFRST",self.mesh_idx_first[idx_surf]+1,slicer=idx_surf) + + # # Reshape the mesh + # # mesh = mesh.ravel(order="C").reshape((3,mesh.shape[0]*mesh.shape[1]), order="F") + # mesh = mesh.transpose((1,0,2)).reshape((mesh.shape[0]*mesh.shape[1],3)) + + # # Set the mesh + # self.set_avl_fort_arr("SURF_MESH_R","MSHBLK",mesh, slicer=(slice(self.mesh_idx_first[idx_surf],self.mesh_idx_first[idx_surf]+nx*ny),slice(0,3))) + + + def set_section_naca(self, isec: int, isurf: int, nasec: int, naca: str, xfminmax: np.ndarray): """Sets the airfoil oml points for the specified surface and section. Computes camber lines, thickness, and oml shape from NACA 4-digit specification. @@ -1210,23 +1385,28 @@ def post_check_input(self, inputDict: dict): # check number of sections, controls, and dvs if len(inputDict["surfaces"]) > 0: surf_names = list(inputDict["surfaces"].keys()) + idx_surf = 0 for i in range(len(inputDict["surfaces"])): surf_dict = inputDict["surfaces"][surf_names[i]] - if self.avl.SURF_GEOM_I.NSEC[i] != surf_dict["num_sections"]: + if self.avl.SURF_GEOM_I.NSEC[idx_surf] != surf_dict["num_sections"]: raise RuntimeError( - f"Mismatch: NSEC[i] = {self.avl.SURF_GEOM_I.NSEC[i]}, Dictionary: {surf_dict['num_sections']}" + f"Mismatch: NSEC[i] = {self.avl.SURF_GEOM_I.NSEC[idx_surf]}, Dictionary: {surf_dict['num_sections']}" ) # Check controls and design variables per section for j in range(surf_dict["num_sections"]): - if self.avl.SURF_GEOM_I.NSCON[i, j] != surf_dict["num_controls"][j]: + if self.avl.SURF_GEOM_I.NSCON[j, idx_surf] != surf_dict["num_controls"][j]: raise RuntimeError( - f"Mismatch: NSCON[i,j] = {self.avl.SURF_GEOM_I.NSCON[i, j]}, Dictionary: {surf_dict['num_controls'][j]}" + f"Mismatch: NSCON[i,j] = {self.avl.SURF_GEOM_I.NSCON[j, idx_surf]}, Dictionary: {surf_dict['num_controls'][j]}" ) - if self.avl.SURF_GEOM_I.NSDES[i, j] != surf_dict["num_design_vars"][j]: + if self.avl.SURF_GEOM_I.NSDES[j, idx_surf] != surf_dict["num_design_vars"][j]: raise RuntimeError( - f"Mismatch: NSDES[i,j] = {self.avl.SURF_GEOM_I.NSDES[i, j]}, Dictionary: {surf_dict['num_design_vars'][j]}" + f"Mismatch: NSDES[i,j] = {self.avl.SURF_GEOM_I.NSDES[j, idx_surf]}, Dictionary: {surf_dict['num_design_vars'][j]}" ) + + idx_surf += 1 + if "yduplicate" in surf_dict: + idx_surf += 1 # Check the global control and design var count if "dname" in inputDict.keys(): @@ -1258,7 +1438,7 @@ def execute_run(self, tol: float = 0.00002): """Run the analysis (equivalent to the AVL command `x` in the OPER menu) Args: - tol: the tolerace of the Newton solver used for timing the aircraft + tol: the tolerace of the Newton solver used for triming the aircraft """ self.set_avl_fort_arr("CASE_R", "EXEC_TOL", tol) self.avl.oper() @@ -3789,6 +3969,234 @@ def plot_geom(self, axes=None): plt.axis("equal") plt.show() + def add_mesh_plot_3d_avl( + self, + axis, + color: str = "black", + mesh_style="--", + mesh_linewidth=0.3, + show_mesh: bool = False, + show_control_points: bool = False + ): + """Adds a 3D plot of the aircraft mesh to a 3D axis + Always plots a flattened mesh from the AVL geometry definition or + a flattened version of any directly assigned meshes. + + Args: + axis: axis to add the plot to + color: what color should the mesh be + mesh_style: line style of the interior mesh, e.g. '-' or '--' + mesh_linewidth: width of the interior mesh, 1.0 will match the surface outline + show_mesh: flag to show the interior mesh of the geometry + """ + mesh_size = self.get_mesh_size() + num_control_surfs = self.get_num_control_surfs() + num_strips = self.get_num_strips() + num_surfs = self.get_num_surfaces() + + # get the mesh points for ploting + mesh_slice = (slice(0, mesh_size),) + strip_slice = (slice(0, num_strips),) + surf_slice = (slice(0, num_surfs),) + + rv1 = self.get_avl_fort_arr("VRTX_R", "RV1", slicer=mesh_slice) # Vortex Left points + rv2 = self.get_avl_fort_arr("VRTX_R", "RV2", slicer=mesh_slice) # Vortex Right points + rc = self.get_avl_fort_arr("VRTX_R", "RC", slicer=mesh_slice) # Control Points + rle1 = self.get_avl_fort_arr("STRP_R", "RLE1", slicer=strip_slice) # Strip left end LE point + rle2 = self.get_avl_fort_arr("STRP_R", "RLE2", slicer=strip_slice) # Strip right end LE point + chord1 = self.get_avl_fort_arr("STRP_R", "CHORD1", slicer=strip_slice) # Left strip chord + chord2 = self.get_avl_fort_arr("STRP_R", "CHORD2", slicer=strip_slice) # Right strip chord + jfrst = self.get_avl_fort_arr("SURF_I", "JFRST", slicer=surf_slice) # Index of first strip in surface + + ijfrst = self.get_avl_fort_arr("STRP_I", "IJFRST", slicer=strip_slice) # Index of first element in strip + nvstrp = self.get_avl_fort_arr("STRP_I", "NVSTRP", slicer=strip_slice) # Number of elements in strip + + nj = self.get_avl_fort_arr("SURF_I", "NJ", slicer=surf_slice) # Number of elements along span in surface + imags = self.get_avl_fort_arr("SURF_I", "IMAGS") # Is surface YDUPL one? + + for idx_surf in range(num_surfs): + # get the range of the elements that belong to this surfaces + strip_st = jfrst[idx_surf] - 1 + strip_end = strip_st + nj[idx_surf] + + # inboard and outboard of outline + # get surfaces that have not been duplicated + if imags[idx_surf] > 0: + j1 = strip_st + jn = strip_end - 1 + dj = 1 + else: + # this surface is a duplicate + j1 = strip_end - 1 + jn = strip_st + dj = -1 + + pts = { + "x": [rle1[j1, 0], rle1[j1, 0] + chord1[j1]], + "y": [rle1[j1, 1], rle1[j1, 1]], + "z": [rle1[j1, 2], rle1[j1, 2]], + } + # # chord-wise grid + axis.plot(pts['x'], pts['y'],pts['z'], color=color) + + pts = { + "x": np.array([rle2[jn, 0], rle2[jn, 0] + chord2[jn]]), + "y": np.array([rle2[jn, 1], rle2[jn, 1]]), + "z": np.array([rle2[jn, 2], rle2[jn, 2]]), + } + + # # chord-wise grid + axis.plot(pts['x'], pts['y'],pts['z'], color=color) + + # # --- outline of surface --- + # front + pts = { + "x": np.append(rle1[j1:jn:dj, 0], rle2[jn, 0]), + "y": np.append(rle1[j1:jn:dj, 1], rle2[jn, 1]), + "z": np.append(rle1[j1:jn:dj, 2], rle2[jn, 2]), + } + axis.plot(pts['x'], pts['y'],pts['z'],"-", color=color) + + # aft + + pts = { + "x": np.append(rle1[j1:jn:dj, 0] + chord1[j1:jn:dj], rle2[jn, 0] + chord2[jn]), + "y": np.append(rle1[j1:jn:dj, 1], rle2[jn, 1]), + "z": np.append(rle1[j1:jn:dj, 2], rle2[jn, 2]), + } + axis.plot(pts['x'], pts['y'],pts['z'],"-", color=color) + + if show_mesh: + for idx_strip in range(strip_st, strip_end): + if ((imags[idx_surf] > 0) and (idx_strip != strip_st)) or ( + (imags[idx_surf] < 0) and (idx_strip != strip_end) + ): + pts = { + "x": [rle1[idx_strip, 0], rle1[idx_strip, 0] + chord1[idx_strip]], + "y": [rle1[idx_strip, 1], rle1[idx_strip, 1]], + "z": [rle1[idx_strip, 2], rle1[idx_strip, 2]], + } + + # # chord-wise grid + axis.plot(pts['x'], pts['y'], pts['z'], mesh_style, color=color, alpha=1.0, linewidth=mesh_linewidth) + + vor_st = ijfrst[idx_strip] - 1 + vor_end = vor_st + nvstrp[idx_strip] + + # spanwise grid + for idx_vor in range(vor_st, vor_end): + pts = { + "x": [rv1[idx_vor, 0], rv2[idx_vor, 0]], + "y": [rv1[idx_vor, 1], rv2[idx_vor, 1]], + "z": [rv1[idx_vor, 2], rv2[idx_vor, 2]], + } + axis.plot(pts['x'], pts['y'],pts['z'], mesh_style, color=color, alpha=1.0, linewidth=mesh_linewidth) + + if show_control_points: + for idx_strip in range(strip_st, strip_end): + vor_st = ijfrst[idx_strip] - 1 + vor_end = vor_st + nvstrp[idx_strip] + + # spanwise grid + for idx_vor in range(vor_st, vor_end): + pts = { + "x": [rc[idx_vor, 0]], + "y": [rc[idx_vor, 1]], + "z": [rc[idx_vor, 2]], + } + axis.scatter(pts['x'], pts['y'], pts['z'],color='r', alpha=1.0, linewidth=0.1) + + def add_mesh_plot_3d_direct( + self, + axis, + color: str = "black", + mesh_style="--", + mesh_linewidth=0.3, + show_mesh: bool = False, + # show_avl_geom: bool = False, + # show_avl_mesh: bool = False, + # avl_mesh_color: str = "red", + # avl_mesh_style: str = "--", + # show_avl_control_points: bool = False + ): + """Plots the true mesh assigned to SURF_MESH AVL common block data on a 3D axis. + Can also plot the mesh stored in SURF on the same axis. + + Args: + axis: axis to add the plot to + color: what color should the mesh be + mesh_style: line style of the interior mesh, e.g. '-' or '--' + mesh_linewidth: width of the interior mesh, 1.0 will match the surface outline + show_mesh: flag to show the interior mesh of the geometry + """ + num_surfs = self.get_num_surfaces() + + for idx_surf in range(num_surfs): + # get the mesh block data + mesh = self.get_mesh(idx_surf) + mesh_x = mesh[:, :, 0] + mesh_y = mesh[:, :, 1] + mesh_z = mesh[:, :, 2] + + # Plot mesh outline + axis.plot(mesh_x[0, :], mesh_y[0, :], mesh_z[0, :], "-", color=color, lw=1) + axis.plot(mesh_x[-1, :], mesh_y[-1, :],mesh_z[-1, :], "-", color=color, lw=1) + axis.plot(mesh_x[:, 0], mesh_y[:, 0],mesh_z[:, 0], "-", color=color, lw=1) + axis.plot(mesh_x[:, -1], mesh_y[:, -1],mesh_z[:, -1], "-", color=color, lw=1) + + if show_mesh: + for i in range(mesh_x.shape[0]): + axis.plot(mesh_x[i, :], mesh_y[i, :], mesh_z[i, :], mesh_style, color=color, lw=mesh_linewidth, alpha=1.0) + for j in range(mesh_x.shape[1]): + axis.plot(mesh_x[:, j], mesh_y[:, j],mesh_z[:, j], mesh_style, color=color, lw=mesh_linewidth, alpha=1.0) + + # if show_avl_geom: + # self.add_mesh_plot_3d_avl( + # axis, + # color = avl_mesh_color, + # mesh_style=avl_mesh_style, + # mesh_linewidth=mesh_linewidth, + # show_mesh = show_avl_mesh, + # show_control_points = show_avl_control_points) + + def plot_geom_3d(self, axes=None, plot_avl_mesh = True, plot_direct_mesh = False): + """Generate a matplotlib plot of geometry + + Args: + axes: Matplotlib axis object to add the plots too. If none are given, the axes will be generated. + """ + + if axes == None: + import matplotlib.pyplot as plt + + ax1 = plt.subplot(projection='3d') + ax1.set_ylabel("X", rotation=0) + ax1.set_aspect("equal") + ax1._axis3don = False + plt.subplots_adjust(left=0.025, right=0.925, top=0.925, bottom=0.025) + else: + ax1, ax2 = axes + + if plot_avl_mesh: + self.add_mesh_plot_3d_avl(ax1, + color = "red", + mesh_style="--", + mesh_linewidth=1.0, + show_mesh= True, + show_control_points = False) + + if plot_direct_mesh: + self.add_mesh_plot_3d_direct(ax1, + color = "black", + mesh_style="--", + mesh_linewidth=0.3, + show_mesh= True) + + if axes == None: + # assume that if we don't provide axes that we want to see the plot + plt.axis("equal") + plt.show() + def get_cp_data(self) -> Tuple[List[np.ndarray], List[np.ndarray]]: """Gets the current surface mesh and cp distribution diff --git a/optvl/utils/check_surface_dict.py b/optvl/utils/check_surface_dict.py index 5784619..cc4c1f2 100755 --- a/optvl/utils/check_surface_dict.py +++ b/optvl/utils/check_surface_dict.py @@ -106,6 +106,10 @@ def pre_check_input_dict(input_dict: dict): "nspans", # number of spanwise elements vector, overriden by nspans "sspaces", # spanwise spacing vector (for each section), overriden by sspace "use surface spacing", # surface spacing set under the surface heeading (known as LSURFSPACING in AVL) + # Geometery: Mesh + "mesh", + "iptloc", + "flatten_mesh", # Control Surfaces # "dname" # IMPLEMENT THIS "icontd", # control variable index @@ -193,7 +197,7 @@ def pre_check_input_dict(input_dict: dict): stacklevel=2, ) input_dict[key] = np.sign(input_dict[key]) - + # Check for keys not implemented if key not in keys_implemented_general: warnings.warn( @@ -201,16 +205,59 @@ def pre_check_input_dict(input_dict: dict): category=RuntimeWarning, stacklevel=2, ) + total_global_control = 0 total_global_design_var = 0 if "surfaces" in input_dict.keys(): if len(input_dict["surfaces"]) > 0: for surface in input_dict["surfaces"].keys(): + # Check if we are directly providing a mesh + if "mesh" in input_dict["surfaces"][surface].keys(): + # Check if sections are specified + if "num_sections" in input_dict["surfaces"][surface].keys(): + # Check if the section indices are provided + if "iptloc" in input_dict["surfaces"][surface].keys(): + # If they are make sure we provide one for every section + if len(input_dict["surfaces"][surface]["iptloc"]) != input_dict["surfaces"][surface]["num_sections"]: + raise ValueError("iptloc vector length does not match num_sections") + # Check if the user provided nspans instead + elif "nspans" in input_dict["surfaces"][surface].keys(): + # setting iptloc to 0 is how we tell the Fortran layer to use nspans + input_dict["surfaces"][surface]["iptloc"] = np.zeros(input_dict["surfaces"][surface]["num_sections"]) + # The OptVL class will have to call the fudging routine to try and auto cut the mesh into sections + else: + warnings.warn( + "Mesh provided for surface dict `{}` for {} sections but locations not defined.\n OptVL will automatically define section locations as close to equally as possible.".format( + surface, input_dict["surfaces"][surface]["num_sections"] + ), + category=RuntimeWarning, + stacklevel=2, + ) + else: + # Assume we have two sections at the ends of mesh and inform the user + warnings.warn( + "Mesh provided for surface dict `{}` but no sections provided.\n Assuming 2 sections at tips.".format( + surface + ), + category=RuntimeWarning, + stacklevel=2, + ) + input_dict["surfaces"][surface]["iptloc"] = np.array([0,input_dict["surfaces"][surface]["mesh"].shape[1]-1],dtype=np.int32) + input_dict["surfaces"][surface]["num_sections"] = 2 + # Verify at least two section if input_dict["surfaces"][surface]["num_sections"] < 2: raise RuntimeError("Must have at least two sections per surface!") + # if no controls are specified then fill it in with 0s + if "num_controls" not in input_dict["surfaces"][surface].keys(): + input_dict["surfaces"][surface]["num_controls"] = np.zeros(input_dict["surfaces"][surface]["num_sections"],dtype=np.int32) + + # if no dvs are specified then fill it in with 0s + if "num_design_vars" not in input_dict["surfaces"][surface].keys(): + input_dict["surfaces"][surface]["num_design_vars"] = np.zeros(input_dict["surfaces"][surface]["num_sections"],dtype=np.int32) + #Checks to see that at most only one of the options in af_load_ops or one of the options in manual_af_override is selected if len(airfoil_spec_keys & input_dict["surfaces"][surface].keys()) > 1: raise RuntimeError( diff --git a/src/aic.f b/src/aic.f index f1453e9..1c8bfae 100644 --- a/src/aic.f +++ b/src/aic.f @@ -119,6 +119,11 @@ SUBROUTINE VVOR(BETM,IYSYM,YSYM,IZSYM,ZSYM,VRCORE, & RV2(1,J),RV2(2,J),RV2(3,J), & BETM,U,V,W,RCORE) C + ! print *, "Influence of", J, "on", I + ! print *, "U:", U + ! print *, "V:", V + ! print *, "W:", W + ! print *, "MARK" IF(IYSYM.NE.0) THEN C... Calculate the influence of the y-IMAGE vortex LBOUND = .TRUE. diff --git a/src/amake.f b/src/amake.f index 8d3fb09..00ea9b2 100644 --- a/src/amake.f +++ b/src/amake.f @@ -186,6 +186,8 @@ SUBROUTINE MAKESURF(ISURF) C C----- fudge spacing array to make nodes match up exactly with interior sections DO ISEC = 2, NSEC(ISURF)-1 + ! Throws an error in the case where the same node is the closest node + ! to two consecutive sections IPT1 = IPTLOC(ISEC-1) IPT2 = IPTLOC(ISEC ) IF(IPT1.EQ.IPT2) THEN @@ -637,10 +639,956 @@ SUBROUTINE MAKESURF(ISURF) C RETURN END ! MAKESURF + + subroutine adjust_mesh_spacing(isurf, nx, ny, mesh, + & iptloc, nsecsurf) + ! This routine is a modified standalone version of the "fudging" + ! operation in makesurf. The main purpose is to deal with cases + ! where the user provide a mesh and does not specify the indicies + ! where the sections are nor do they include the number of spanwise + ! elements associated with each section. This routine is intended + ! to be run as a preprocessing step to compute iptloc and the fudged mesh + ! as once we have iptloc makesurf_mesh will know how to handle the sections. + INCLUDE 'AVL.INC' + ! input/output + integer nx, ny, nsecsurf, isurf + integer isec, ipt, ipt1, ipt2, idx_sec, idx_pt + integer iptloc(nsecsurf) + real mesh(3,nx,ny) + + ! working variables + integer niptloc + real ylen(nsecsurf), yzlen(nsecsurf) + real yptloc, yptdel, yp1, yp2, dy, dz, y_mesh, dy_mesh + + ! check that iptloc is the correct size + if (nsecsurf /= NSEC(isurf)) then + write(*,'(A,I2,A,I2)') 'given size of iptloc:',nsecsurf, + & ' does not match NSEC(isurf):', NSEC(isurf) + endif + + + ! Check if the mesh can be adjusted + if (ny < nsecsurf) then + print *, "*** Not enought spanwise nodes to split the mesh" + stop + end if + + ! Unlike the standard fudging routine we have no idea where + ! each section's leading edge is located ahead of time + ! Instead we have to make an initial guess by cutting up + ! the mesh into equal pieces spanwise assuming the wing is flat. + ! We only need to do this is there isn't already a guess for iptloc + + if (iptloc(1) .eq. 0) then + ! compute mesh y length + y_mesh = mesh(2,1,ny) - mesh(2,1,1) + dy_mesh = y_mesh/(NSEC(isurf)-1) + write(*,*) 'y_mesh', y_mesh + write(*,*) 'dy_mesh', dy_mesh + + ! Chop up into equal y length pieces + ylen(1) = 0. + do idx_sec = 2,NSEC(isurf) + ylen(idx_sec) = ylen(idx_sec-1) + dy_mesh + end do + + ! Find node nearest each section + do idx_sec = 2, NSEC(isurf)-1 + yptloc = 1.0E9 + iptloc(idx_sec) = 1 + do idx_pt = 1, ny + yptdel = abs(mesh(2,1,1)+ ylen(idx_sec) - mesh(2,1,idx_pt)) + if(yptdel .LT. yptloc) then + yptloc = yptdel + iptloc(idx_sec) = idx_pt + endif + enddo + enddo + iptloc(1) = 1 + iptloc(NSEC(isurf)) = ny + end if + + ! NOTE-SB: I don't think we need this + ! I originally included it to be more consistent with AVL + ! The prior routine only considers the y distance while + ! AVL considers y and z. However, the above routine appears + ! to work fine on its own and running this after appears to + ! cause issues. + + ! Now compute yz arc length using the computed section indicies +! yzlen(1) = 0. +! do idx_sec = 2, NSEC(isurf) +! dy = mesh(2,1,iptloc(idx_sec)) - mesh(2,1 +! & , (iptloc(idx_sec)-1)) +! dz = mesh(3,1,iptloc(idx_sec)) - mesh(3,1 +! & , (iptloc(idx_sec)-1)) +! yzlen(idx_sec) = yzlen(idx_sec-1) + sqrt(dy*dy + dz*dz) +! end do + +! ! Now do the AVL fudging routine to ensure the sections don't split panels + +! ! Find node nearest each section +! do isec = 2, NSEC(isurf)-1 +! yptloc = 1.0E9 +! iptloc(isec) = 1 +! do ipt = 1, ny +! yptdel = abs(yzlen(isec) - mesh(2,1,ipt)) +! if(yptdel .LT. yptloc) then +! yptloc = yptdel +! iptloc(ISEC) = ipt +! endif +! enddo +! enddo +! iptloc(1) = 1 +! iptloc(NSEC(ISURF)) = ny + +! print *, "Final iptloc", iptloc + + ! fudge spacing array to make nodes match up exactly with interior sections + do isec = 2, NSEC(isurf)-1 + ! Throws an error in the case where the same node is the closest node + ! to two consecutive sections + ipt1 = iptloc(isec-1) + ipt2 = iptloc(isec ) + if(ipt1.EQ.ipt2) then + CALL STRIP(STITLE(isurf),NST) + WRITE(*,7000) isec, STITLE(isurf)(1:NST) + STOP + end if + + ! fudge spacing to this section so that nodes match up exactly with section + ypt1 = mesh(2,1,ipt1) + yscale = (yzlen(isec)-yzlen(isec-1)) / (mesh(2,1,ipt2) + & -ypt1) + do ipt = ipt1, ipt2-1 + mesh(2,1,ipt) = yzlen(isec-1) + yscale*(mesh(2,1,ipt)-ypt1) + end do + + ! check for unique spacing node for next section, if not we need more nodes + ipt1 = iptloc(isec ) + ipt2 = iptloc(isec+1) + if(ipt1.EQ.ipt2) then + CALL STRIP(STITLE(isurf),NST) + WRITE(*,7000) isec, STITLE(isurf)(1:NST) + STOP + endif + + ! fudge spacing to this section so that nodes match up exactly with section + ypt1 = mesh(2,1,ipt1) + ypt2 = mesh(2,1,ipt2) + yscale = (ypt2-yzlen(isec)) / (ypt2-ypt1) + do ipt = ipt1, ipt2-1 + mesh(2,1,ipt) = yzlen(isec) + yscale*(mesh(2,1,ipt)-ypt1) + enddo + + 7000 format( + & /' *** Cannot adjust spanwise spacing at section', I3, + & ', on surface ', A + & /' *** Insufficient number of spanwise vortices to work with') + enddo + + end subroutine adjust_mesh_spacing + + integer function flatidx(idx_x, idx_y, idx_surf) + include 'AVL.INC' + ! store MFRST and NVC in the common block + integer idx_x, idx_y, idx_surf + flatidx = idx_x + (idx_y - 1) * (NVC(idx_surf)+1) + return + end function flatidx + + subroutine makesurf_mesh(isurf) +c-------------------------------------------------------------- +c Sets up all stuff for surface ISURF, +C using info from configuration input +C and the given mesh coordinate array. +c-------------------------------------------------------------- + INCLUDE 'AVL.INC' + ! input/output + integer isurf + + ! working variables (AVL original) + PARAMETER (KCMAX=50, + & KSMAX=500) + REAL CHSIN, CHCOS, CHSINL, CHSINR, CHCOSL, CHCOSR, AINCL, AINCR, + & CHORDL, CHORDR, CLAFL, CLAFR, SLOPEL, SLOPER, DXDX, ZU_L, + & ZL_L, ZU_R, ZL_R, ZL, ZR, SUM, WTOT, ASTRP + REAL CHSINL_G(NGMAX),CHCOSL_G(NGMAX), + & CHSINR_G(NGMAX),CHCOSR_G(NGMAX), + & XLED(NDMAX), XTED(NDMAX), GAINDA(NDMAX) + INTEGER ISCONL(NDMAX), ISCONR(NDMAX) + + ! working variables (OptVL additions) + real m1, m2, m3, f1, f2, fc, dc1, dc2, dc, a1, a2, a3, xptxind1, + & xptxind2 + real mesh_surf(3,(NVC(isurf)+1)*(NVS(isurf)+1)) + integer idx_vor, idx_strip, idx_sec, idx_dim, idx_coef, idx_x, + & idx_node, idx_nodel, idx_noder, idx_node_yp1, idx_node_nx, + & idx_node_nx_yp1, idx_y, nx, ny + + ! functions + integer flatidx + + ! Get data from common block + nx = NVC(isurf) + 1 + ny = NVS(isurf) + 1 + + ! If the user doesn't input a index vector telling us at what + ! spanwise index each section is located they will have to have + ! provided nspans otherwise they will have to go back and provide + ! iptloc or run adjust_mesh_spacing as a preprocessing step to get + ! a iptloc vector. + if (IPTSEC(1,isurf) .eq. 0) then + ! if NSPANS is given then use it + if (NSPANS(1,isurf) .ne. 0) then + IPTSEC(1,isurf) = 1 + do idx_sec = 2,NSEC(isurf) + IPTSEC(idx_sec,isurf) = IPTSEC(idx_sec-1,isurf) + + & NSPANS(idx_sec-1,isurf) + end do + else + print *, '* Provide NSPANS or IPTSEC. (Hint: Run adjust_mesh_& + & spacing)' + stop + end if + end if + + ! Check MFRST + if (MFRST(isurf) .eq. 0) then + print *, "* Provide the index where the mesh begins for surface", + & isurf + end if + + ! Get the mesh for this surface from the the common block + mesh_surf = MSHBLK(:,MFRST(isurf):MFRST(isurf)+(nx*ny)-1) + + ! Perform input checks from makesurf + + IF(NSEC(ISURF).LT.2) THEN + WRITE(*,*) '*** Need at least 2 sections per surface.' + STOP + ENDIF + + IF(NVC(ISURF).GT.KCMAX) THEN + WRITE(*,*) '* makesurf_mesh: Array overflow. Increase KCMAX to', + & NVC(ISURF) + NVC(ISURF) = KCMAX + ENDIF + + IF(NVS(ISURF).GT.KSMAX) THEN + WRITE(*,*) '* makesurf_mesh: Array overflow. Increase KSMAX to', + & NVS(ISURF) + NVS(ISURF) = KSMAX + ENDIF + + ! Image flag set to indicate section definition direction + ! IMAGS= 1 defines edge 1 located at surface root edge + ! IMAGS=-1 defines edge 2 located at surface root edge (reflected surfaces) + IMAGS(ISURF) = 1 + + ! Start accumulating the element and strip index references + ! Accumulate the first element in surface + if (ISURF == 1) then + IFRST(ISURF) = 1 + else + IFRST(ISURF) = IFRST(ISURF-1) + NK(ISURF-1)*NJ(ISURF-1) + endif + + ! Accumulate the first strip in surface + if (ISURF == 1) then + JFRST(ISURF) = 1 + else + JFRST(ISURF) = JFRST(ISURF-1) + NJ(ISURF-1) + endif + + ! Set NK from input data (python layer will ensure this is consistent) + NK(ISURF) = NVC(ISURF) + + ! We need to start counting strips now since it's a global count + idx_strip = JFRST(ISURF) + + ! Bypass the entire spanwise node generation routine and go straight to store counters + ! Index of first section in surface + IF (ISURF .EQ. 1) THEN + ICNTFRST(ISURF) = 1 + ELSE + ICNTFRST(ISURF) = ICNTFRST(ISURF-1) + NCNTSEC(ISURF-1) + ENDIF + ! Number of sections in surface + NCNTSEC(ISURF) = NSEC(ISURF) + ! Store the spanwise index of each section in each surface + DO ISEC = 1, NSEC(ISURF) + II = ICNTFRST(ISURF) + (ISEC-1) + ICNTSEC(II) = IPTSEC(ISEC,isurf) + ENDDO + + + ! Apply the scaling and translations to the mesh as a whole + do idx_y = 1,ny + do idx_x = 1,nx + do idx_dim = 1,3 + idx_node = flatidx(idx_x, idx_y, isurf) + mesh_surf(idx_dim,idx_node) = XYZSCAL(idx_dim,isurf) + & *mesh_surf(idx_dim,idx_node) + XYZTRAN(idx_dim,isurf) + end do + end do + end do + + + ! Setup the strips + + ! Set spanwise elements to 0 + NJ(ISURF) = 0 + + ! Check control and design vars (input routine should've already checked this tbh) + IF(NCONTROL.GT.NDMAX) THEN + WRITE(*,*) '*** Too many control variables. Increase NDMAX to', + & NCONTROL + STOP + ENDIF + + IF(NDESIGN.GT.NGMAX) THEN + WRITE(*,*) '*** Too many design variables. Increase NGMAX to', + & NDESIGN + STOP + ENDIF + + + ! Loop over sections + do idx_sec = 1, NSEC(isurf)-1 + + ! Set reference information for the section + iptl = IPTSEC(idx_sec,isurf) + iptr = IPTSEC(idx_sec+1,isurf) + nspan = iptr - iptl + NJ(isurf) = NJ(isurf) + nspan + + + ! We need to compute the chord and claf values at the left and right edge of the section + ! These will be needed by AVL for control surface setup and control point placement + idx_node = flatidx(1,iptl,isurf) + idx_node_nx = flatidx(nx,iptl,isurf) + CHORDL = sqrt((mesh_surf(1,idx_node_nx)-mesh_surf(1,idx_node))**2 + & + (mesh_surf(3,idx_node_nx)-mesh_surf(3,idx_node))**2) + idx_node = flatidx(1,iptr,isurf) + idx_node_nx = flatidx(nx,iptr,isurf) + CHORDR = sqrt((mesh_surf(1,idx_node_nx)-mesh_surf(1,idx_node))**2 + & + (mesh_surf(3,idx_node_nx)-mesh_surf(3,idx_node))**2) + CLAFL = CLAF(idx_sec, isurf) + CLAFR = CLAF(idx_sec+1,isurf) + + ! Compute the incidence angle at the section end points + ! We will need this later to iterpolate chord projections + ! SAB Note: This type of interpolation assumes the section + ! is linear. However, the twist angles it produces can be applied + ! to an arbitrary mesh. The user just needs to be aware of what they are + ! applying here. + ! Analogy: imagine that we create a trapazoidal section + ! that matches up with the root and tip chords of your arbitrary wing section + ! The trapzoid can encompass the wing section or parts of the section can be + ! protruding out side the trapazoid. It doesn't matter. Now we twist the trapazoid + ! so that the angles at the root and tip match what is specified at the sections in + ! AINCS. However, when we twist we make sure to keep the leading and trailing edges + ! linear (straight line along the LE and TE). The angles at each strip required to + ! do are what gets applied to the normal vector at each strip. + AINCL = AINCS(idx_sec,isurf)*DTR + ADDINC(isurf)*DTR + AINCR = AINCS(idx_sec+1,isurf)*DTR + ADDINC(isurf)*DTR + CHSINL = CHORDL*SIN(AINCL) + CHSINR = CHORDR*SIN(AINCR) + CHCOSL = CHORDL*COS(AINCL) + CHCOSR = CHORDR*COS(AINCR) + + ! We need to determine which controls belong to this section + ! Bring over the routine for this from makesurf + DO N = 1, NCONTROL + ISCONL(N) = 0 + ISCONR(N) = 0 + DO ISCON = 1, NSCON(idx_sec,isurf) + IF(ICONTD(ISCON,idx_sec,isurf) .EQ.N) ISCONL(N) = ISCON + ENDDO + DO ISCON = 1, NSCON(idx_sec+1,isurf) + IF(ICONTD(ISCON,idx_sec+1,isurf).EQ.N) ISCONR(N) = ISCON + ENDDO + ENDDO + + ! We need to determine which dvs belong to this section + ! and setup the chord projection gains + ! Bring over the routine for this from makesurf + DO N = 1, NDESIGN + CHSINL_G(N) = 0. + CHSINR_G(N) = 0. + CHCOSL_G(N) = 0. + CHCOSR_G(N) = 0. + + DO ISDES = 1, NSDES(idx_sec,isurf) + IF(IDESTD(ISDES,idx_sec,isurf).EQ.N) THEN + CHSINL_G(N) = CHCOSL * GAING(ISDES,idx_sec,isurf)*DTR + CHCOSL_G(N) = -CHSINL * GAING(ISDES,idx_sec,isurf)*DTR + ENDIF + ENDDO + + DO ISDES = 1, NSDES(idx_sec+1,isurf) + IF(IDESTD(ISDES,idx_sec+1,isurf).EQ.N) THEN + CHSINR_G(N) = CHCOSR * GAING(ISDES,idx_sec+1,isurf)*DTR + CHCOSR_G(N) = -CHSINR * GAING(ISDES,idx_sec+1,isurf)*DTR + ENDIF + ENDDO + ENDDO + + + ! Set the strip geometry data + ! Note these computations assume the mesh is not necessarily planar + ! ultimately if/when we flatten the mesh into a planar one we will want + ! to use the leading edge positions and chords from the original input mesh + + ! Loop over strips in section + do ispan = 1,nspan + idx_y = idx_strip - JFRST(isurf) + 1 + + ! Strip left side + idx_node = flatidx(1,idx_y,isurf) + idx_node_nx = flatidx(nx,idx_y,isurf) + do idx_dim = 1,3 + RLE1(idx_dim,idx_strip) = mesh_surf(idx_dim,idx_node) + end do + CHORD1(idx_strip) = sqrt((mesh_surf(1,idx_node_nx) + & -mesh_surf(1,idx_node))**2 + (mesh_surf(3,idx_node_nx) + & -mesh_surf(3,idx_node))**2) + + ! Strip right side + idx_node_yp1 = flatidx(1,idx_y+1,isurf) + idx_node_nx_yp1 = flatidx(nx,idx_y+1,isurf) + do idx_dim = 1,3 + RLE2(idx_dim,idx_strip) = mesh_surf(idx_dim,idx_node_yp1) + end do + CHORD2(idx_strip) = sqrt((mesh_surf(1,idx_node_nx_yp1) + & -mesh_surf(1,idx_node_yp1))**2 + (mesh_surf(3,idx_node_nx_yp1) + & -mesh_surf(3,idx_node_yp1))**2) + + ! Strip mid-point + do idx_dim = 1,3 + ! Since the strips are linear SPANWISE we can just interpolate + RLE(idx_dim,idx_strip) = (RLE1(idx_dim,idx_strip) + & + RLE2(idx_dim,idx_strip))/2. + end do + ! The strips are not necessarily linear chord wise but by definition the chord value is + ! so we can interpolate + CHORD(idx_strip) = (CHORD1(idx_strip)+CHORD2(idx_strip))/2. + + ! Strip geometric incidence angle at the mid-point + ! This is strip incidence angle is computed from the LE and TE points + ! of the given geometry and is completely independent of AINC + ! This quantity is needed to correctly handle nonplanar meshes and is only needed if the mesh isn't flattened + GINCSTRIP(idx_strip) = atan2(((mesh_surf(3,idx_node_nx) + & + mesh_surf(3,idx_node_nx_yp1))/2.- (mesh_surf(3,idx_node) + + & mesh_surf(3,idx_node_yp1))/2.), + & ((mesh_surf(1,idx_node_nx) + mesh_surf(1,idx_node_nx_yp1))/2. + & - (mesh_surf(1,idx_node) + mesh_surf(1,idx_node_yp1))/2.)) + + ! Strip width + m2 = mesh_surf(2,idx_node_yp1)-mesh_surf(2,idx_node) + m3 = mesh_surf(3,idx_node_yp1)-mesh_surf(3,idx_node) + WSTRIP(idx_strip) = sqrt(m2**2 + m3**2) + + ! Strip LE and TE sweep slopes + tanle(idx_strip) = (mesh_surf(1,idx_node_yp1) + & -mesh_surf(1,idx_node))/WSTRIP(idx_strip) + idx_node = flatidx(nx,idx_y,isurf) + idx_node_yp1 = flatidx(nx,idx_y+1,isurf) + tante(idx_strip) = (mesh_surf(1,idx_node_yp1) + & -mesh_surf(1,idx_node))/WSTRIP(idx_strip) + + ! Compute chord projections and strip twists + ! In AVL the AINCS are not interpolated. The chord projections are + ! So we have to replicate this effect. + + ! LINEAR interpolation over the section: left, right, and midpoint + idx_nodel = flatidx(1,iptl,isurf) + idx_noder = flatidx(1,iptr,isurf) + + f1 = (mesh_surf(2,idx_node)-mesh_surf(2,idx_nodel))/ + & (mesh_surf(2,idx_noder)-mesh_surf(2,idx_nodel)) + f2 = (mesh_surf(2,idx_node_yp1)-mesh_surf(2,idx_nodel))/ + & (mesh_surf(2,idx_noder)-mesh_surf(2,idx_nodel)) + fc = (((mesh_surf(2,idx_node_yp1)+mesh_surf(2,idx_node))/2.) + & -mesh_surf(2,idx_nodel))/(mesh_surf(2,idx_noder) + & -mesh_surf(2,idx_nodel)) + + ! Strip left side incidence + CHSIN = CHSINL + f1*(CHSINR-CHSINL) + CHCOS = CHCOSL + f1*(CHCOSR-CHCOSL) + AINC1(idx_strip) = ATAN2(CHSIN,CHCOS) + + ! Strip right side incidence + CHSIN = CHSINL + f2*(CHSINR-CHSINL) + CHCOS = CHCOSL + f2*(CHCOSR-CHCOSL) + AINC2(idx_strip) = ATAN2(CHSIN,CHCOS) + + ! Strip mid-point incidence + CHSIN = CHSINL + fc*(CHSINR-CHSINL) + CHCOS = CHCOSL + fc*(CHCOSR-CHCOSL) + AINC(idx_strip) = ATAN2(CHSIN,CHCOS) + + ! Set dv gains for incidence angles + ! Bring over the routine for this from make surf + DO N = 1, NDESIGN + CHSIN_G = (1.0-FC)*CHSINL_G(N) + FC*CHSINR_G(N) + CHCOS_G = (1.0-FC)*CHCOSL_G(N) + FC*CHCOSR_G(N) + AINC_G(idx_strip,N) = (CHCOS*CHSIN_G - CHSIN*CHCOS_G) + & / (CHSIN**2 + CHCOS**2) + ENDDO + + ! We have to now setup any control surfaces we defined for this section + ! Bring over the routine for this from makesurf + DO N = 1, NCONTROL + ICL = ISCONL(N) + ICR = ISCONR(N) + + IF(ICL.EQ.0 .OR. ICR.EQ.0) THEN + ! no control effect here + GAINDA(N) = 0. + XLED(N) = 0. + XTED(N) = 0. + + VHINGE(1,idx_strip,N) = 0. + VHINGE(2,idx_strip,N) = 0. + VHINGE(3,idx_strip,N) = 0. + + VREFL(idx_strip,N) = 0. + + PHINGE(1,idx_strip,N) = 0. + PHINGE(2,idx_strip,N) = 0. + PHINGE(3,idx_strip,N) = 0. + + ELSE + ! control variable # N is active here + GAINDA(N) = GAIND(ICL,idx_sec ,isurf)*(1.0-FC) + & + GAIND(ICR,idx_sec+1,isurf)* FC + + ! SAB Note: This interpolation ensures that the hinge line is + ! is linear which I think it is an ok assumption for arbitrary wings as long as the user is aware + ! A curve hinge line could work if needed if we just interpolate XHINGED and scaled by local chord + XHD = CHORDL*XHINGED(ICL,idx_sec ,isurf)*(1.0-FC) + & + CHORDR*XHINGED(ICR,idx_sec+1,isurf)* FC + IF(XHD.GE.0.0) THEN + ! TE control surface, with hinge at XHD + XLED(N) = XHD + XTED(N) = CHORD(idx_strip) + ELSE + ! LE control surface, with hinge at -XHD + XLED(N) = 0.0 + XTED(N) = -XHD + ENDIF + + VHX = VHINGED(1,ICL,idx_sec,isurf)*XYZSCAL(1,isurf) + VHY = VHINGED(2,ICL,idx_sec,isurf)*XYZSCAL(2,isurf) + VHZ = VHINGED(3,ICL,idx_sec,isurf)*XYZSCAL(3,isurf) + VSQ = VHX**2 + VHY**2 + VHZ**2 + IF(VSQ.EQ.0.0) THEN + ! default: set hinge vector along hingeline + ! We are just setting the hinge line across the section + ! this assumes the hinge is linear even for a nonlinear wing + VHX = mesh_surf(1,idx_noder) + & + ABS(CHORDR*XHINGED(ICR,idx_sec+1,isurf)) + & - mesh_surf(1,idx_nodel) + & - ABS(CHORDL*XHINGED(ICL,idx_sec,isurf)) + VHY = mesh_surf(2,idx_noder) + & - mesh_surf(2,idx_nodel) + VHZ = mesh_surf(3,idx_noder) + & - mesh_surf(3,idx_nodel) + VHX = VHX*XYZSCAL(1,isurf) + VHY = VHY*XYZSCAL(2,isurf) + VHZ = VHZ*XYZSCAL(3,isurf) + VSQ = VHX**2 + VHY**2 + VHZ**2 + ENDIF + + VMOD = SQRT(VSQ) + VHINGE(1,idx_strip,N) = VHX/VMOD + VHINGE(2,idx_strip,N) = VHY/VMOD + VHINGE(3,idx_strip,N) = VHZ/VMOD + + VREFL(idx_strip,N) = REFLD(ICL,idx_sec, isurf) + + IF(XHD .GE. 0.0) THEN + PHINGE(1,idx_strip,N) = RLE(1,idx_strip) + XHD + PHINGE(2,idx_strip,N) = RLE(2,idx_strip) + PHINGE(3,idx_strip,N) = RLE(3,idx_strip) + ELSE + PHINGE(1,idx_strip,N) = RLE(1,idx_strip) - XHD + PHINGE(2,idx_strip,N) = RLE(2,idx_strip) + PHINGE(3,idx_strip,N) = RLE(3,idx_strip) + ENDIF + ENDIF + ENDDO + + ! Interpolate CD-CL polar defining data from input sections to strips + DO idx_coef = 1, 6 + CLCD(idx_coef,idx_strip) = (1.0-fc)* + & CLCDSEC(idx_coef,idx_sec,isurf) + + & fc*CLCDSEC(idx_coef,idx_sec+1,isurf) + END DO + ! If the min drag is zero flag the strip as no-viscous data + LVISCSTRP(idx_strip) = (CLCD(4,idx_strip).NE.0.0) + + + ! Set the panel (vortex) geometry data + + ! Accumulate the strip element indicies and start counting vorticies + if (idx_strip .eq. 1) then + IJFRST(idx_strip) = 1 + else + IJFRST(idx_strip) = IJFRST(idx_strip - 1) + + & NVSTRP(idx_strip - 1) + endif + idx_vor = IJFRST(idx_strip) + NVSTRP(idx_strip) = NVC(isurf) + + ! Associate the strip with the surface + NSURFS(idx_strip) = isurf + + ! Prepare for cross section interpolation + NSL = NASEC(idx_sec , isurf) + NSR = NASEC(idx_sec+1, isurf) + + ! CHORDC = CHORD(idx_strip) + + ! Interpolate claf over the section + ! SAB: In AVL this quantity is interpolated as a product with chord + ! We then divide by the chord at the strip to recover claf at the strip + ! This only works correctly for linear sections. For arbitrary sections + ! this can result in the claf varying across the span even when the claf + ! between two secions is equal. + ! After reaching out to Hal Youngren it is determined that it is + ! best to just interpolate claf straight up for now +! clafc = (1.-FC)*(CHORDL/CHORD(idx_strip))*CLAFL +! & + FC *(CHORDR/CHORD(idx_strip))*CLAFR + clafc = (1.-fc)*clafl + fc*clafr + + ! loop over vorticies for the strip + do idx_x = 1, nvc(isurf) + + ! Left bound vortex points + idx_node = flatidx(idx_x,idx_y,isurf) + ! Compute the panel's left side chord + dc1 = sqrt((mesh_surf(1,idx_node+1) - mesh_surf(1,idx_node))**2 + & + (mesh_surf(3,idx_node+1) - mesh_surf(3,idx_node))**2) + + if (LMESHFLAT(isurf)) then + ! Place vortex at panel quarter chord of the flat mesh + dx1 = sqrt((mesh_surf(1,idx_node) - RLE1(1,idx_strip))**2 + & + (mesh_surf(3,idx_node) - RLE1(3,idx_strip))**2) + RV1(2,idx_vor) = RLE1(2,idx_strip) + RV1(3,idx_vor) = RLE1(3,idx_strip) + RV1(1,idx_vor) = RLE1(1,idx_strip) + dx1 + (dc1/4.) + + ! Compute the panel's left side angle + a1 = atan2((mesh_surf(3,idx_node+1) - mesh_surf(3,idx_node)), + & (mesh_surf(1,idx_node+1) - mesh_surf(1,idx_node))) + ! Place vortex at panel quarter chord of the true mesh + RV1MSH(2,idx_vor) = mesh_surf(2,idx_node) + RV1MSH(1,idx_vor) = mesh_surf(1,idx_node) + (dc1/4.)*cos(a1) + RV1MSH(3,idx_vor) = mesh_surf(3,idx_node) + (dc1/4.)*sin(a1) + else + ! Compute the panel's left side angle + a1 = atan2((mesh_surf(3,idx_node+1) - mesh_surf(3,idx_node)), + & (mesh_surf(1,idx_node+1) - mesh_surf(1,idx_node))) + ! Place vortex at panel quarter chord + RV1(2,idx_vor) = mesh_surf(2,idx_node) + RV1(1,idx_vor) = mesh_surf(1,idx_node) + (dc1/4.)*cos(a1) + RV1(3,idx_vor) = mesh_surf(3,idx_node) + (dc1/4.)*sin(a1) + + ! Make a copy in the true mesh array for post processing + RV1MSH(2,idx_vor) = RV1(2,idx_vor) + RV1MSH(1,idx_vor) = RV1(1,idx_vor) + RV1MSH(3,idx_vor) = RV1(3,idx_vor) + end if + + ! Right bound vortex points + idx_node_yp1 = flatidx(idx_x,idx_y+1,isurf) + ! Compute the panel's right side chord + dc2 = sqrt((mesh_surf(1,idx_node_yp1+1) + & - mesh_surf(1,idx_node_yp1))**2 + (mesh_surf(3,idx_node_yp1+1) + & - mesh_surf(3,idx_node_yp1))**2) + + if (LMESHFLAT(isurf)) then + ! Place vortex at panel quarter chord of the flat mesh + dx2 = sqrt((mesh_surf(1,idx_node_yp1) - RLE2(1,idx_strip))**2 + & + (mesh_surf(3,idx_node_yp1) - RLE2(3,idx_strip))**2) + RV2(2,idx_vor) = RLE2(2,idx_strip) + RV2(3,idx_vor) = RLE2(3,idx_strip) + RV2(1,idx_vor) = RLE2(1,idx_strip) + dx2 + (dc2/4.) + + ! Compute the panel's right side angle + a2 = atan2((mesh_surf(3,idx_node_yp1+1) - + & mesh_surf(3,idx_node_yp1)), (mesh_surf(1,idx_node_yp1+1) - + & mesh_surf(1,idx_node_yp1))) + ! Place vortex at panel quarter chord of the true mesh + RV2MSH(2,idx_vor) = mesh_surf(2,idx_node_yp1) + RV2MSH(1,idx_vor) = mesh_surf(1,idx_node_yp1) + (dc2/4.)*cos(a2) + RV2MSH(3,idx_vor) = mesh_surf(3,idx_node_yp1) + (dc2/4.)*sin(a2) + else + ! Compute the panel's right side angle + a2 = atan2((mesh_surf(3,idx_node_yp1+1) - + & mesh_surf(3,idx_node_yp1)), (mesh_surf(1,idx_node_yp1+1) - + & mesh_surf(1,idx_node_yp1))) + ! Place vortex at panel quarter chord + RV2(2,idx_vor) = mesh_surf(2,idx_node_yp1) + RV2(1,idx_vor) = mesh_surf(1,idx_node_yp1) + (dc2/4.)*cos(a2) + RV2(3,idx_vor) = mesh_surf(3,idx_node_yp1) + (dc2/4.)*sin(a2) + + ! Make a copy in the true mesh array for post processing + RV2MSH(2,idx_vor) = RV2(2,idx_vor) + RV2MSH(1,idx_vor) = RV2(1,idx_vor) + RV2MSH(3,idx_vor) = RV2(3,idx_vor) + end if + + ! Mid-point bound vortex points + ! Compute the panel's mid-point chord + ! Panels themselves can never be curved so just interpolate the chord + ! store as the panel chord in common block + DXV(idx_vor) = (dc1+dc2)/2. + ! We need to compute the midpoint angle and panel strip chord projection + ! as we need them to compute normals based on the real mesh + a3 = atan2(((mesh_surf(3,idx_node_yp1+1) + & + mesh_surf(3,idx_node+1))/2.- (mesh_surf(3,idx_node_yp1) + + & mesh_surf(3,idx_node))/2.), + & ((mesh_surf(1,idx_node_yp1+1) + mesh_surf(1,idx_node+1))/2. + & - (mesh_surf(1,idx_node_yp1) + mesh_surf(1,idx_node))/2.)) + ! project the panel chord onto the strip chord + DXSTRPV(idx_vor) = DXV(idx_vor)*cos(a3-GINCSTRIP(idx_strip)) + + if (LMESHFLAT(isurf)) then + ! Place vortex at panel quarter chord of the flat mesh + dx3 = sqrt(((mesh_surf(1,idx_node_yp1)+mesh_surf(1,idx_node))/2 + & - RLE(1,idx_strip))**2 + & + ((mesh_surf(3,idx_node_yp1)+mesh_surf(3,idx_node))/2 + & - RLE(3,idx_strip))**2) + RV(2,idx_vor) = RLE(2,idx_strip) + RV(3,idx_vor) = RLE(3,idx_strip) + RV(1,idx_vor) = RLE(1,idx_strip) + dx3 + (DXV(idx_vor)/4.) + + ! Place vortex at panel quarter chord of the true mesh + RVMSH(2,idx_vor) = (mesh_surf(2,idx_node_yp1) + & + mesh_surf(2,idx_node))/2. + RVMSH(1,idx_vor) = (mesh_surf(1,idx_node_yp1) + & +mesh_surf(1,idx_node))/2.+ (DXV(idx_vor)/4.)*cos(a3) + RVMSH(3,idx_vor) = (mesh_surf(3,idx_node_yp1) + & +mesh_surf(3,idx_node))/2. + (DXV(idx_vor)/4.)*sin(a3) + else + ! Place vortex at panel quarter chord + RV(2,idx_vor) = (mesh_surf(2,idx_node_yp1) + & + mesh_surf(2,idx_node))/2. + RV(1,idx_vor) = (mesh_surf(1,idx_node_yp1) + & +mesh_surf(1,idx_node))/2.+ (DXV(idx_vor)/4.)*cos(a3) + RV(3,idx_vor) = (mesh_surf(3,idx_node_yp1) + & +mesh_surf(3,idx_node))/2. + (DXV(idx_vor)/4.)*sin(a3) + + ! Make a copy in the true mesh array for post processing + RVMSH(2,idx_vor) = RV(2,idx_vor) + RVMSH(1,idx_vor) = RV(1,idx_vor) + RVMSH(3,idx_vor) = RV(3,idx_vor) + end if + + + ! Panel Control points + ! Y- point + ! is just the panel midpoint + RC(2,idx_vor) = RV(2,idx_vor) + ! Place the control point at the quarter chord + half chord*clafc + ! note that clafc is a scaler so is 1. is for 2pi + ! use data from vortex mid-point computation + if (LMESHFLAT(isurf)) then + RC(1,idx_vor) = RV(1,idx_vor) + clafc*(DXV(idx_vor)/2.) + RC(3,idx_vor) = RV(3,idx_vor) + + RCMSH(1,idx_vor) = RVMSH(1,idx_vor) + & + clafc*(DXV(idx_vor)/2.)*cos(a3) + RCMSH(3,idx_vor) = RVMSH(3,idx_vor) + & + clafc*(DXV(idx_vor)/2.)*sin(a3) + RCMSH(2,idx_vor) = RVMSH(2,idx_vor) + else + RC(1,idx_vor) = RV(1,idx_vor) + clafc*(DXV(idx_vor)/2.)*cos(a3) + RC(3,idx_vor) = RV(3,idx_vor) + clafc*(DXV(idx_vor)/2.)*sin(a3) + + ! Make a copy in the true mesh array for post processing + RCMSH(1,idx_vor) = RC(1,idx_vor) + RCMSH(3,idx_vor) = RC(3,idx_vor) + RCMSH(2,idx_vor) = RC(2,idx_vor) + end if + + ! Source points + ! Y- point + RS(2,idx_vor) = RV(2,idx_vor) + ! Place the source point at the half chord + ! use data from vortex mid-point computation + ! add another quarter chord to the quarter chord + if (LMESHFLAT(isurf)) then + RS(1,idx_vor) = RV(1,idx_vor) + (DXV(idx_vor)/4.) + RS(3,idx_vor) = RV(3,idx_vor) + (DXV(idx_vor)/4.) + else + RS(1,idx_vor) = RV(1,idx_vor) + (DXV(idx_vor)/4.)*cos(a3) + RS(3,idx_vor) = RV(3,idx_vor) + (DXV(idx_vor)/4.)*sin(a3) + end if + + + ! Set the camber slopes for the panel + + ! Camber slope at control point + CALL AKIMA(XASEC(1,idx_sec, isurf),SASEC(1,idx_sec, isurf), + & NSL,(RC(1,idx_vor)-RLE(1,idx_strip)) + & /CHORD(idx_strip),SLOPEL, DSDX) + CALL AKIMA(XASEC(1,idx_sec+1,isurf),SASEC(1,idx_sec+1,isurf), + & NSR,(RC(1,idx_vor)-RLE(1,idx_strip)) + & /CHORD(idx_strip),SLOPER, DSDX) + + ! Interpolate this as is per Hal Youngren (for now) + SLOPEC(idx_vor) = (1.-fc)*SLOPEL + fc*SLOPER +! SLOPEC(idx_vor) = (1.-fc)*(CHORDL/CHORD(idx_strip))*SLOPEL +! & + fc *(CHORDR/CHORD(idx_strip))*SLOPER + + ! Camber slope at vortex mid-point + CALL AKIMA(XASEC(1,idx_sec, isurf),SASEC(1,idx_sec, isurf), + & NSL,(RV(1,idx_vor)-RLE(1,idx_strip)) + & /CHORD(idx_strip),SLOPEL, DSDX) + CALL AKIMA(XASEC(1,idx_sec+1,isurf),SASEC(1,idx_sec+1,isurf), + & NSR,(RV(1,idx_vor)-RLE(1,idx_strip)) + & /CHORD(idx_strip),SLOPER, DSDX) + + ! Interpolate this as is per Hal Youngren (for now) + SLOPEV(idx_vor) = (1.-fc)*SLOPEL + fc*SLOPER +! SLOPEV(idx_vor) = (1.-fc)*(CHORDL/CHORD(idx_strip))*SLOPEL +! & + fc *(CHORDR/CHORD(idx_strip))*SLOPER + + ! Associate the panel with it's strip's chord and component + CHORDV(idx_vor) = CHORD(idx_strip) + NSURFV(idx_vor) = LSCOMP(isurf) + + ! Enforce no penetration at the control point + LVNC(idx_vor) = .true. + + ! element inherits alpha,beta flag from surface + LVALBE(idx_vor) = LFALBE(isurf) + + ! We need to scale the control surface gains by the fraction + ! of the element on the control surface + do N = 1, NCONTROL + !scale control gain by factor 0..1, (fraction of element on control surface) + FRACLE = (XLED(N)/CHORD(idx_strip)-((mesh_surf(1,idx_node_yp1) + & -mesh_surf(1,idx_node))/2.)/CHORD(idx_strip)) / + & (DXV(idx_vor)/CHORD(idx_strip)) + + FRACTE = (XTED(N)/CHORD(idx_strip)-((mesh_surf(1,idx_node_yp1) + & -mesh_surf(1,idx_node))/2.)/CHORD(idx_strip)) / + & (DXV(idx_vor)/CHORD(idx_strip)) + + FRACLE = MIN( 1.0 , MAX( 0.0 , FRACLE ) ) + FRACTE = MIN( 1.0 , MAX( 0.0 , FRACTE ) ) + + DCONTROL(idx_vor,N) = GAINDA(N)*(FRACTE-FRACLE) + end do + + ! TE control point used only if surface sheds a wake + LVNC(idx_vor) = LFWAKE(isurf) + + ! Use the cross sections to generate the OML + ! nodal grid associated with vortex strip (aft-panel nodes) + ! NOTE: airfoil in plane of wing, but not rotated perpendicular to dihedral; + ! retained in (x,z) plane at this point + + ! Store the panel LE mid point for the next panel in the strip + ! This gets used a lot here + ! We use the original input mesh (true mesh) to compute points for the OML + xptxind1 = ((mesh_surf(1,idx_node+1)+mesh_surf(1,idx_node_yp1+1)) + & /2 - RLE(1,idx_strip))/CHORD(idx_strip) + +! xptxind2 = (mesh_surf(1,idx_node_yp1+1) +! & - RLE2(1,idx_strip))/CHORD2(idx_strip) + + ! Interpolate cross section on left side + CALL AKIMA( XLASEC(1,idx_sec,isurf), ZLASEC(1,idx_sec,isurf), + & NSL,xptxind1, ZL_L, DSDX ) + CALL AKIMA( XUASEC(1,idx_sec,isurf), ZUASEC(1,idx_sec,isurf), + & NSL,xptxind1, ZU_L, DSDX ) + + ! Interpolate cross section on right side + CALL AKIMA( XLASEC(1,idx_sec+1,isurf),ZLASEC(1,idx_sec+1,isurf), + & NSR, xptxind1, ZL_R, DSDX) + + CALL AKIMA( XUASEC(1,idx_sec+1,isurf),ZUASEC(1,idx_sec+1,isurf), + & NSR, xptxind1, ZU_R, DSDX) + + + ! Compute the left aft node of panel + ! X-point + XYN1(1,idx_vor) = RLE1(1,idx_strip) + + & xptxind1*CHORD1(idx_strip) + + ! Y-point + XYN1(2,idx_vor) = RLE1(2,idx_strip) + + ! Interpolate z from sections to left aft node of panel + ZL = (1.-f1)*ZL_L + f1 *ZL_R + ZU = (1.-f1)*ZU_L + f1 *ZU_R + + ! Store left aft z-point + ZLON1(idx_vor) = RLE1(3,idx_strip) + ZL*CHORD1(idx_strip) + ZUPN1(idx_vor) = RLE1(3,idx_strip) + ZU*CHORD1(idx_strip) + + ! Compute the right aft node of panel + ! X-point + XYN2(1,idx_vor) = RLE2(1,idx_strip) + + & xptxind1*CHORD2(idx_strip) + + ! Y-point + XYN2(2,idx_vor) = RLE2(2,idx_strip) + + ! Interpolate z from sections to right aft node of panel + ZL = (1.-f2)*ZL_L + f2 *ZL_R + ZU = (1.-f2)*ZU_L + f2 *ZU_R + + ! Store right aft z-point + ZLON2(idx_vor) = RLE2(3,idx_strip) + ZL*CHORD2(idx_strip) + ZUPN2(idx_vor) = RLE2(3,idx_strip) + ZU*CHORD2(idx_strip) + + + idx_vor = idx_vor + 1 + end do ! End vortex loop + idx_strip = idx_strip + 1 + end do ! End strip loop + + end do ! End section loop + + ! Compute the wetted area and cave from the true mesh + sum = 0.0 + wtot = 0.0 + DO JJ = 1, NJ(isurf) + J = JFRST(isurf) + JJ-1 + ASTRP = WSTRIP(J)*CHORD(J) + SUM = SUM + ASTRP + WTOT = WTOT + WSTRIP(J) + ENDDO + SSURF(isurf) = SUM + + IF(WTOT .EQ. 0.0) THEN + CAVESURF(isurf) = 0.0 + ELSE + CAVESURF(isurf) = sum/wtot + ENDIF + ! add number of strips to the global count + NSTRIP = NSTRIP + NJ(isurf) + ! add number of of votrices to the global count + NVOR = NVOR + NK(isurf)*NJ(isurf) + + end subroutine makesurf_mesh subroutine update_surfaces() c-------------------------------------------------------------- c Updates all surfaces, using the stored data. +c Resets the strips and vorticies so that AVL rebuilds them +c from the updated geometry. Recomputes panels normals and +c tells AVL to rebuild the AIC and other aero related data +c arrays on the next execution. c-------------------------------------------------------------- include 'AVL.INC' @@ -659,9 +1607,17 @@ subroutine update_surfaces() ! it was probably duplicated from the previous one cycle end if - call makesurf(ISURF) + if (lsurfmsh(isurf)) then + call makesurf_mesh(ISURF) + else + call makesurf(ISURF) + end if else - call makesurf(ISURF) + if (lsurfmsh(isurf)) then + call makesurf_mesh(ISURF) + else + call makesurf(ISURF) + end if endif if(ldupl(isurf)) then @@ -671,11 +1627,11 @@ subroutine update_surfaces() CALL ENCALC - LAIC = .FALSE. - LSRD = .FALSE. - LVEL = .FALSE. - LSOL = .FALSE. - LSEN = .FALSE. + LAIC = .FALSE. ! Tell AVL that the AIC is no longer valid and to regenerate it + LSRD = .FALSE. ! Tell AVL that unit source+doublet strengths are no longer valid and to regenerate them + LVEL = .FALSE. ! Tell AVL that the induced velocity matrix is no longer valid and to regenerate it + LSOL = .FALSE. ! Tell AVL that a valid solution no longer exists + LSEN = .FALSE. ! Tell AVL that valid sensitives no longer exist end subroutine update_surfaces @@ -941,6 +1897,8 @@ SUBROUTINE SDUPL(NN, Ypt,MSG) LFLOAD(NNI) = LFLOAD(NN) LRANGE(NNI) = LRANGE(NN) LSURFSPACING(NNI) = LSURFSPACING(NN) + LMESHFLAT(NNI) = LMESHFLAT(NN) + LSURFMSH(NNI) = LSURFMSH(NN) C---- accumulate stuff for new image surface ! IFRST(NNI) = NVOR + 1 @@ -1000,6 +1958,7 @@ SUBROUTINE SDUPL(NN, Ypt,MSG) RLE(2,JJI) = -RLE(2,JJ) + YOFF RLE(3,JJI) = RLE(3,JJ) CHORD(JJI) = CHORD(JJ) + GINCSTRIP(JJI) = GINCSTRIP(JJ) WSTRIP(JJI) = WSTRIP(JJ) TANLE(JJI) = -TANLE(JJ) AINC (JJI) = AINC(JJ) @@ -1065,10 +2024,27 @@ SUBROUTINE SDUPL(NN, Ypt,MSG) SLOPEC(III) = SLOPEC(II) SLOPEV(III) = SLOPEV(II) DXV(III) = DXV(II) + DXSTRPV(III) = DXSTRPV(II) CHORDV(III) = CHORDV(II) NSURFV(III) = LSCOMP(NNI) LVALBE(III) = LVALBE(II) LVNC(III) = LVNC(II) + ! Duplicate mesh data if we are using a mesh + if (lsurfmsh(NN)) then + RV1MSH(1,III) = RV2MSH(1,II) + RV1MSH(2,III) = -RV2MSH(2,II) + YOFF + RV1MSH(3,III) = RV2MSH(3,II) + RV2MSH(1,III) = RV1MSH(1,II) + RV2MSH(2,III) = -RV1MSH(2,II) + YOFF + RV2MSH(3,III) = RV1MSH(3,II) + RVMSH(1,III) = RVMSH(1,II) + RVMSH(2,III) = -RVMSH(2,II) + YOFF + RVMSH(3,III) = RVMSH(3,II) + RCMSH(1,III) = RCMSH(1,II) + RCMSH(2,III) = -RCMSH(2,II) + YOFF + RCMSH(3,III) = RCMSH(3,II) + end if + C DO N = 1, NCONTROL ccc RSGN = SIGN( 1.0 , VREFL(JJ,N) ) @@ -1176,6 +2152,8 @@ SUBROUTINE ENCALC C...PURPOSE To calculate the normal vectors for the strips, C the horseshoe vortices, and the control points. C Incorporates surface deflections. + ! Also checks if surface has been assigned a point cloud mesh + ! and uses the real mesh to compute normals if it is C C...INPUT NVOR Number of vortices C X1 Coordinates of endpoint #1 of the vortices @@ -1199,10 +2177,46 @@ SUBROUTINE ENCALC C REAL EP(3), EQ(3), ES(3), EB(3), EC(3), ECXB(3) REAL EC_G(3,NDMAX), ECXB_G(3) + real(kind=avl_real) :: dchstrip, DXT, DYT, DZT C C...Calculate the normal vector at control points and bound vortex midpoints C DO 10 J = 1, NSTRIP + + ! Since we cannot seperate the encalc routine for direct mesh assignment we have to make it a branch here + if (lsurfmsh(nsurfs(J))) then + + ! Calculate normal vector for the strip (normal to X axis) + ! we can't just interpolate this anymore given that + ! the strip is no longer necessarily linear chordwise + + ! We want the spanwise unit vector for the strip at the + ! chordwise location specified by SAXFR (usually set to 0.25) + ! Loop over all panels in the strip until we find the one that contains + ! the SAXFR position in it's projected chord. Since the panels themselves are still linear + ! we can just use the bound vortex unit vector of that panel as + ! the spanwise unit vector of the strip at SAXFR + + ! SAB: This is slow, find a better way to do this + dchstrip = 0.0 + searchSAXFR: do i = IJFRST(J),IJFRST(J) + (NVSTRP(J)-1) + dchstrip = dchstrip+DXSTRPV(i) + if (dchstrip .ge. CHORD(J)*SAXFR) then + exit searchSAXFR + end if + end do searchSAXFR + + + ! compute the spanwise unit vector for Vperp def + DXT = RV2MSH(1,I)-RV1MSH(1,I) + DYT = RV2MSH(2,I)-RV1MSH(2,I) + DZT = RV2MSH(3,I)-RV1MSH(3,I) + XSREF(J) = RVMSH(1,I) + YSREF(J) = RVMSH(2,I) + ZSREF(J) = RVMSH(3,I) + + else + ! original encalc routine for standard AVL geometry C C...Calculate normal vector for the strip (normal to X axis) I = IJFRST(J) @@ -1230,19 +2244,21 @@ SUBROUTINE ENCALC DXT = (1.0-SAXFR)*DXLE + SAXFR*DXTE DYT = (1.0-SAXFR)*DYLE + SAXFR*DYTE DZT = (1.0-SAXFR)*DZLE + SAXFR*DZTE -C - ESS(1,J) = DXT/SQRT(DXT*DXT + DYT*DYT + DZT*DZT) - ESS(2,J) = DYT/SQRT(DXT*DXT + DYT*DYT + DZT*DZT) - ESS(3,J) = DZT/SQRT(DXT*DXT + DYT*DYT + DZT*DZT) -C - ENSY(J) = -DZT/SQRT(DYT*DYT + DZT*DZT) - ENSZ(J) = DYT/SQRT(DYT*DYT + DZT*DZT) C XSREF(J) = (1.0-SAXFR)*AXLE + SAXFR*AXTE YSREF(J) = (1.0-SAXFR)*AYLE + SAXFR*AYTE ZSREF(J) = (1.0-SAXFR)*AZLE + SAXFR*AZTE + end if C C + ESS(1,J) = DXT/SQRT(DXT*DXT + DYT*DYT + DZT*DZT) + ESS(2,J) = DYT/SQRT(DXT*DXT + DYT*DYT + DZT*DZT) + ESS(3,J) = DZT/SQRT(DXT*DXT + DYT*DYT + DZT*DZT) + + ! Treffz plane normals + ENSY(J) = -DZT/SQRT(DYT*DYT + DZT*DZT) + ENSZ(J) = DYT/SQRT(DYT*DYT + DZT*DZT) + ES(1) = 0. ES(2) = ENSY(J) ES(3) = ENSZ(J) @@ -1271,11 +2287,19 @@ SUBROUTINE ENCALC ENC_G(2,I,N) = 0. ENC_G(3,I,N) = 0. ENDDO + + if (lsurfmsh(nsurfs(J))) then + ! Define unit vector along bound leg + DXB = RV2MSH(1,I)-RV1MSH(1,I) ! right h.v. pt - left h.v. pt + DYB = RV2MSH(2,I)-RV1MSH(2,I) + DZB = RV2MSH(3,I)-RV1MSH(3,I) + else C C...Define unit vector along bound leg DXB = RV2(1,I)-RV1(1,I) ! right h.v. pt - left h.v. pt DYB = RV2(2,I)-RV1(2,I) DZB = RV2(3,I)-RV1(3,I) + end if EMAG = SQRT(DXB**2 + DYB**2 + DZB**2) EB(1) = DXB/EMAG EB(2) = DYB/EMAG @@ -1284,6 +2308,11 @@ SUBROUTINE ENCALC C...Define direction of normal vector at control point C The YZ projection of the normal vector matches the camber slope C + section local incidence in the YZ defining plane for the section + ! First start by combining the contributions to the panel + ! incidence from AVL incidence and camberline slope variables + ! these are not actual geometric transformations of the mesh + ! but rather further modifications to the chordwise vector that + ! will get used to compute normals ANG = AINC(J) - ATAN(SLOPEC(I)) cc IF(LDES) THEN C--------- add design-variable contribution to angle @@ -1294,10 +2323,32 @@ SUBROUTINE ENCALC C SINC = SIN(ANG) COSC = COS(ANG) + + if (lsurfmsh(nsurfs(J))) then + ! direct mesh assignemnt branch + ! now we compute the chordwise panel vector + ! note that panel's chordwise vector has contributions + ! from both the geometry itself and the incidence modification + ! from the AVL AINC and camber slope variables + + ! To avoid storing uncessary info in the common block + ! Get the geometric chordwise vector using RVMSH and RCMSH which should + ! be located in the same plane given that each individual panel is a + ! plane + + ! Note that like in AVL the sin of the incidence is projected + ! to the strip's normal in the YZ plane (Treffz plane) + ! which is ES(2) and ES(3) computed earlier + EC(1) = COSC + (RCMSH(1,I)-RVMSH(1,I)) + EC(2) = -SINC*ES(2) + (RCMSH(2,I)-RVMSH(2,I)) + EC(3) = -SINC*ES(3) + (RCMSH(3,I)-RVMSH(3,I)) + + else EC(1) = COSC EC(2) = -SINC*ES(2) EC(3) = -SINC*ES(3) ! EC = rotation of strip normal vector? or along chord? + end if DO N = 1, NDESIGN EC_G(1,N) = -SINC *AINC_G(J,N) EC_G(2,N) = -COSC*ES(2)*AINC_G(J,N) @@ -1330,6 +2381,11 @@ SUBROUTINE ENCALC C...Define direction of normal vector at vortex mid-point. C The YZ projection of the normal vector matches the camber slope C + section local incidence in the YZ defining plane for the section + ! This section is identical to the normal vector at the control + ! point. The only different is that the AVL camberline slope + ! is taken at the bound vortex point rather than the control point + ! the geometric contributions to the normal vector at both of these + ! point is identical as the lie in the plane of the same panel. ANG = AINC(J) - ATAN(SLOPEV(I)) cc IF(LDES) THEN C--------- add design-variable contribution to angle @@ -1340,9 +2396,17 @@ SUBROUTINE ENCALC C SINC = SIN(ANG) COSC = COS(ANG) + if (lsurfmsh(nsurfs(J))) then + ! direct mesh assignment branch + EC(1) = COSC + (RCMSH(1,I)-RVMSH(1,I)) + EC(2) = -SINC*ES(2) + (RCMSH(2,I)-RVMSH(2,I)) + EC(3) = -SINC*ES(3) + (RCMSH(3,I)-RVMSH(3,I)) + else EC(1) = COSC EC(2) = -SINC*ES(2) EC(3) = -SINC*ES(3) + end if + DO N = 1, NDESIGN EC_G(1,N) = -SINC *AINC_G(J,N) EC_G(2,N) = -COSC*ES(2)*AINC_G(J,N) @@ -1376,6 +2440,8 @@ SUBROUTINE ENCALC C C======================================================= C-------- rotate normal vectors for control surface + ! this is a pure rotation of the normal vector + ! the geometric contribution from the mesh is already accounted for DO 100 N = 1, NCONTROL C C---------- skip everything if this element is unaffected by control variable N @@ -1438,3 +2504,303 @@ SUBROUTINE ENCALC RETURN END ! ENCALC + + + + +! SUBROUTINE ENCALCMSH +! C +! C...PURPOSE To calculate the normal vectors for the strips, +! C the horseshoe vortices, and the control points. +! C Assuming arbitrary point cloud mesh +! C Incorporates surface deflections. +! C +! C...INPUT NVOR Number of vortices +! C X1 Coordinates of endpoint #1 of the vortices +! C X2 Coordinates of endpoint #2 of the vortices +! C SLOPEV Slope at bound vortices +! C SLOPEC Slope at control points +! C NSTRIP Number of strips +! C IJFRST Index of first element in strip +! C NVSTRP No. of vortices in strip +! C AINC Angle of incidence of strip +! C LDES include design-variable deflections if TRUE +! C +! C...OUTPUT ENC(3) Normal vector at control point +! C ENV(3) Normal vector at bound vortices +! C ENSY, ENSZ Strip normal vector (ENSX=0) +! C LSTRIPOFF Non-used strip (T) (below z=ZSYM) +! C +! C...COMMENTS +! C +! INCLUDE 'AVL.INC' +! C +! REAL EP(3), EQ(3), ES(3), EB(3), EC(3), ECXB(3) +! REAL EC_G(3,NDMAX), ECXB_G(3) + +! real(kind=avl_real) :: dchstrip, DXT, DYT, DZT +! C +! C...Calculate the normal vector at control points and bound vortex midpoints +! C +! DO 10 J = 1, NSTRIP +! C +! C...Calculate normal vector for the strip (normal to X axis) +! ! we can't just interpolate this anymore given that +! ! the strip is no longer necessarily linear chordwise + +! ! We want the spanwise unit vector for the strip at the +! ! chordwise location specified by SAXFR (usually set to 0.25) +! ! Loop over all panels in the strip until we find the one that contains +! ! the SAXFR position in it's projected chord. Since the panels themselves are still linear +! ! we can just use the bound vortex unit vector of that panel as +! ! the spanwise unit vector of the strip at SAXFR + +! ! SAB: This is slow, find a better way to do this +! dchstrip = 0.0 +! searchSAXFR: do i = IJFRST(J),IJFRST(J) + (NVSTRP(J)-1) +! dchstrip = dchstrip+DXSTRPV(i) +! if (dchstrip .ge. CHORD(J)*SAXFR) then +! exit searchSAXFR +! end if +! end do searchSAXFR + +! ! print *, "I", I + +! ! compute the spanwise unit vector for Vperp def +! DXT = RV2MSH(1,I)-RV1MSH(1,I) +! DYT = RV2MSH(2,I)-RV1MSH(2,I) +! DZT = RV2MSH(3,I)-RV1MSH(3,I) +! XSREF(J) = RVMSH(1,I) +! YSREF(J) = RVMSH(2,I) +! ZSREF(J) = RVMSH(3,I) + +! ! print *, "DVT", DYT +! ! print *, "RV2(2,I)-RV1(2,I)", RV2(2,I)-RV1(2,I) +! ! print *, "RV2(2,I)", RV2(2,I) +! ! print *, "RV1(2,I)", RV1(2,I) +! ! print *, "NSTRIP", NSTRIP +! ! print *, "J", J +! ESS(1,J) = DXT/SQRT(DXT*DXT + DYT*DYT + DZT*DZT) +! ESS(2,J) = DYT/SQRT(DXT*DXT + DYT*DYT + DZT*DZT) +! ESS(3,J) = DZT/SQRT(DXT*DXT + DYT*DYT + DZT*DZT) + +! ! Treffz plane normals +! ENSY(J) = -DZT/SQRT(DYT*DYT + DZT*DZT) +! ENSZ(J) = DYT/SQRT(DYT*DYT + DZT*DZT) + +! ES(1) = 0. +! ES(2) = ENSY(J) +! ES(3) = ENSZ(J) +! C +! LSTRIPOFF(J) = .FALSE. +! C +! NV = NVSTRP(J) +! DO 105 II = 1, NV +! C +! I = IJFRST(J) + (II-1) +! C +! DO N = 1, NCONTROL +! ENV_D(1,I,N) = 0. +! ENV_D(2,I,N) = 0. +! ENV_D(3,I,N) = 0. +! ENC_D(1,I,N) = 0. +! ENC_D(2,I,N) = 0. +! ENC_D(3,I,N) = 0. +! ENDDO +! C +! DO N = 1, NDESIGN +! ENV_G(1,I,N) = 0. +! ENV_G(2,I,N) = 0. +! ENV_G(3,I,N) = 0. +! ENC_G(1,I,N) = 0. +! ENC_G(2,I,N) = 0. +! ENC_G(3,I,N) = 0. +! ENDDO +! C +! C...Define unit vector along bound leg +! DXB = RV2MSH(1,I)-RV1MSH(1,I) ! right h.v. pt - left h.v. pt +! DYB = RV2MSH(2,I)-RV1MSH(2,I) +! DZB = RV2MSH(3,I)-RV1MSH(3,I) +! EMAG = SQRT(DXB**2 + DYB**2 + DZB**2) +! EB(1) = DXB/EMAG +! EB(2) = DYB/EMAG +! EB(3) = DZB/EMAG +! C +! C...Define direction of normal vector at control point + +! ! First start by combining the contributions to the panel +! ! incidence from AVL incidence and camberline slope variables +! ! these are not actual geometric transformations of the mesh +! ! but rather further modifications to the chordwise vector that +! ! will get used to compute normals +! ANG = AINC(J) - ATAN(SLOPEC(I)) +! C--------- add design-variable contribution to angle +! DO N = 1, NDESIGN +! ANG = ANG + AINC_G(J,N)*DELDES(N) +! ENDDO +! C +! ! now we compute the chordwise panel vector +! ! note that panel's chordwise vector has contributions +! ! from both the geometry itself and the incidence modification +! ! from the AVL AINC and camber slope variables + +! ! To avoid storing uncessary info in the common block +! ! Get the geometric chordwise vector using RV and RC which should +! ! be located in the same plane given that each individual panel is a +! ! plane + +! ! Note that like in AVL the sin of the incidence is projected +! ! to the strip's normal in the YZ plane (Treffz plane) +! ! which is ES(2) and ES(3) computed earlier +! SINC = SIN(ANG) +! COSC = COS(ANG) +! EC(1) = COSC + (RCMSH(1,I)-RVMSH(1,I)) +! EC(2) = -SINC*ES(2) + (RCMSH(2,I)-RVMSH(2,I)) +! EC(3) = -SINC*ES(3) + (RCMSH(3,I)-RVMSH(3,I)) + +! DO N = 1, NDESIGN +! EC_G(1,N) = -SINC *AINC_G(J,N) +! EC_G(2,N) = -COSC*ES(2)*AINC_G(J,N) +! EC_G(3,N) = -COSC*ES(3)*AINC_G(J,N) +! ENDDO +! C +! C...Normal vector is perpendicular to camberline vector and to the bound leg +! CALL CROSS(EC,EB,ECXB) +! EMAG = SQRT(ECXB(1)**2 + ECXB(2)**2 + ECXB(3)**2) +! IF(EMAG.NE.0.0) THEN +! ENC(1,I) = ECXB(1)/EMAG +! ENC(2,I) = ECXB(2)/EMAG +! ENC(3,I) = ECXB(3)/EMAG +! DO N = 1, NDESIGN +! CALL CROSS(EC_G(1,N),EB,ECXB_G) +! EMAG_G = ENC(1,I)*ECXB_G(1) +! & + ENC(2,I)*ECXB_G(2) +! & + ENC(3,I)*ECXB_G(3) +! ENC_G(1,I,N) = (ECXB_G(1) - ENC(1,I)*EMAG_G)/EMAG +! ENC_G(2,I,N) = (ECXB_G(2) - ENC(2,I)*EMAG_G)/EMAG +! ENC_G(3,I,N) = (ECXB_G(3) - ENC(3,I)*EMAG_G)/EMAG +! ENDDO +! ELSE +! ENC(1,I) = ES(1) +! ENC(2,I) = ES(2) +! ENC(3,I) = ES(3) +! ENDIF + +! C +! C +! C...Define direction of normal vector at vortex mid-point. + +! ! This section is identical to the normal vector at the control +! ! point. The only different is that the AVL camberline slope +! ! is taken at the bound vortex point rather than the control point +! ! the geometric contributions to the normal vector at both of these +! ! point is identical as the lie in the plane of the same panel. +! ANG = AINC(J) - ATAN(SLOPEV(I)) + +! C--------- add design-variable contribution to angle +! DO N = 1, NDESIGN +! ANG = ANG + AINC_G(J,N)*DELDES(N) +! ENDDO + +! C +! SINC = SIN(ANG) +! COSC = COS(ANG) +! EC(1) = COSC + (RCMSH(1,I)-RVMSH(1,I)) +! EC(2) = -SINC*ES(2) + (RCMSH(2,I)-RVMSH(2,I)) +! EC(3) = -SINC*ES(3) + (RCMSH(3,I)-RVMSH(3,I)) +! DO N = 1, NDESIGN +! EC_G(1,N) = -SINC *AINC_G(J,N) +! EC_G(2,N) = -COSC*ES(2)*AINC_G(J,N) +! EC_G(3,N) = -COSC*ES(3)*AINC_G(J,N) +! ENDDO +! C +! C...Normal vector is perpendicular to camberline vector and to the bound leg +! CALL CROSS(EC,EB,ECXB) +! EMAG = SQRT(ECXB(1)**2 + ECXB(2)**2 + ECXB(3)**2) +! IF(EMAG.NE.0.0) THEN +! ENV(1,I) = ECXB(1)/EMAG +! ENV(2,I) = ECXB(2)/EMAG +! ENV(3,I) = ECXB(3)/EMAG +! DO N = 1, NDESIGN +! CALL CROSS(EC_G(1,N),EB,ECXB_G) +! EMAG_G = ENC(1,I)*ECXB_G(1) +! & + ENC(2,I)*ECXB_G(2) +! & + ENC(3,I)*ECXB_G(3) +! ENV_G(1,I,N) = (ECXB_G(1) - ENV(1,I)*EMAG_G)/EMAG +! ENV_G(2,I,N) = (ECXB_G(2) - ENV(2,I)*EMAG_G)/EMAG +! ENV_G(3,I,N) = (ECXB_G(3) - ENV(3,I)*EMAG_G)/EMAG +! ENDDO +! ELSE +! ENV(1,I) = ES(1) +! ENV(2,I) = ES(2) +! ENV(3,I) = ES(3) +! ENDIF +! C +! C +! ccc write(*,*) i, dcontrol(i,1), dcontrol(i,2) +! C +! C======================================================= +! C-------- rotate normal vectors for control surface +! ! this is a pure rotation of the normal vector +! ! the geometric contribution from the mesh is already accounted for +! DO 100 N = 1, NCONTROL +! C +! C---------- skip everything if this element is unaffected by control variable N +! IF(DCONTROL(I,N).EQ.0.0) GO TO 100 +! C +! ANG = DTR*DCONTROL(I,N)*DELCON(N) +! ANG_DDC = DTR*DCONTROL(I,N) +! C +! COSD = COS(ANG) +! SIND = SIN(ANG) +! C +! C---------- EP = normal-vector component perpendicular to hinge line +! ENDOT = DOT(ENC(1,I),VHINGE(1,J,N)) +! EP(1) = ENC(1,I) - ENDOT*VHINGE(1,J,N) +! EP(2) = ENC(2,I) - ENDOT*VHINGE(2,J,N) +! EP(3) = ENC(3,I) - ENDOT*VHINGE(3,J,N) +! C---------- EQ = unit vector perpendicular to both EP and hinge line +! CALL CROSS(VHINGE(1,J,N),EP,EQ) +! C +! C---------- rotated vector would consist of sin,cos parts from EP and EQ, +! C- with hinge-parallel component ENDOT restored +! cc ENC(1,I) = EP(1)*COSD + EQ(1)*SIND + ENDOT*VHINGE(1,J,N) +! cc ENC(2,I) = EP(2)*COSD + EQ(2)*SIND + ENDOT*VHINGE(2,J,N) +! cc ENC(3,I) = EP(3)*COSD + EQ(3)*SIND + ENDOT*VHINGE(3,J,N) +! C +! C---------- linearize about zero deflection (COSD=1, SIND=0) +! ENC_D(1,I,N) = ENC_D(1,I,N) + EQ(1)*ANG_DDC +! ENC_D(2,I,N) = ENC_D(2,I,N) + EQ(2)*ANG_DDC +! ENC_D(3,I,N) = ENC_D(3,I,N) + EQ(3)*ANG_DDC +! C +! C +! C---------- repeat for ENV vector +! C +! C---------- EP = normal-vector component perpendicular to hinge line +! ENDOT = DOT(ENV(1,I),VHINGE(1,J,N)) +! EP(1) = ENV(1,I) - ENDOT*VHINGE(1,J,N) +! EP(2) = ENV(2,I) - ENDOT*VHINGE(2,J,N) +! EP(3) = ENV(3,I) - ENDOT*VHINGE(3,J,N) +! C---------- EQ = unit vector perpendicular to both EP and hinge line +! CALL CROSS(VHINGE(1,J,N),EP,EQ) +! C +! C---------- rotated vector would consist of sin,cos parts from EP and EQ, +! C- with hinge-parallel component ENDOT restored +! cc ENV(1,I) = EP(1)*COSD + EQ(1)*SIND + ENDOT*VHINGE(1,J,N) +! cc ENV(2,I) = EP(2)*COSD + EQ(2)*SIND + ENDOT*VHINGE(2,J,N) +! cc ENV(3,I) = EP(3)*COSD + EQ(3)*SIND + ENDOT*VHINGE(3,J,N) +! C +! C---------- linearize about zero deflection (COSD=1, SIND=0) +! ENV_D(1,I,N) = ENV_D(1,I,N) + EQ(1)*ANG_DDC +! ENV_D(2,I,N) = ENV_D(2,I,N) + EQ(2)*ANG_DDC +! ENV_D(3,I,N) = ENV_D(3,I,N) + EQ(3)*ANG_DDC +! 100 CONTINUE +! 101 CONTINUE +! C +! 105 CONTINUE +! 10 CONTINUE +! C +! LENC = .TRUE. +! C +! RETURN +! END ! ENCALCMSH \ No newline at end of file diff --git a/src/avl.f b/src/avl.f index 5f6e413..0d2b8ca 100644 --- a/src/avl.f +++ b/src/avl.f @@ -460,11 +460,11 @@ SUBROUTINE loadGEO(geom_file) C----- initialize state C CALL VARINI C - LAIC = .FALSE. - LSRD = .FALSE. - LVEL = .FALSE. - LSOL = .FALSE. - LSEN = .FALSE. + LAIC = .FALSE. ! Tell AVL that the AIC is no longer valid and to regenerate it + LSRD = .FALSE. ! Tell AVL that unit source+doublet strengths are no longer valid and to regenerate them + LVEL = .FALSE. ! Tell AVL that the induced velocity matrix is no longer valid and to regenerate it + LSOL = .FALSE. ! Tell AVL that a valid solution no longer exists + LSEN = .FALSE. ! Tell AVL that valid sensitives no longer exists END diff --git a/src/f2py/libavl.pyf b/src/f2py/libavl.pyf index 8960ffb..50c9845 100644 --- a/src/f2py/libavl.pyf +++ b/src/f2py/libavl.pyf @@ -135,6 +135,16 @@ python module libavl ! in real*8, intent(inout) :: asys(jemax,jemax) end subroutine get_system_matrix + subroutine makesurf_mesh(isurf) ! in :libavl:amake.f + integer :: isurf + end subroutine makesurf_mesh + + subroutine adjust_mesh_spacing(isurf, nx, ny, mesh, iptloc, nsecsurf) ! in :libavl:amake.f + integer :: isurf, nx, ny, nsecsurf + integer, intent(inout) :: iptloc(nsecsurf) + real*8 :: mesh(3,nx,ny) + end subroutine adjust_mesh_spacing + !subroutine getcam(x, y, n, xc, yc, tc, nc, lnorm) ! in :libavl:airutil.f !integer, intent(in) :: n !integer, intent(inout) :: nc @@ -168,5 +178,12 @@ python module libavl ! in real*8 :: xb(nb), yb(nb) logical :: storecoords end subroutine set_body_coordinates + + !subroutine update_surface_mesh_HACK(isurf,mesh,nx,ny,iptloc,nsecsurf,lcount,lcall) ! in :libavl:amake.f + !integer :: isurf, nx, ny, nsecsurf + !logical :: lcount, lcall + !real*8 :: mesh(3,nx,ny) + !integer :: iptloc(nsecsurf) + !end subroutine update_surfaces_mesh end interface end python module libavl diff --git a/src/includes/AVL.INC.in b/src/includes/AVL.INC.in index 57a389f..e436046 100644 --- a/src/includes/AVL.INC.in +++ b/src/includes/AVL.INC.in @@ -497,8 +497,30 @@ c & REFLD(ICONX, NSMAX, NFMAX), ! control surface reflection & GAING(ICONX, NSMAX, NFMax) ! desgin variable gain + COMMON /SURF_MESH_I/ + & MFRST(NFMAX), ! stores the index in the MSHBLK where each surface's mesh begins + & IPTSEC(NSMAX,NFMAX) ! stores the iptloc vector for each surface -C !!--- end added variables for python geometry minipulation --- + REAL(kind=avl_real) MSHBLK + REAL(kind=avl_real) RV1MSH + REAL(kind=avl_real) RV2MSH + REAL(kind=avl_real) RVMSH + REAL(kind=avl_real) RCMSH + COMMON /SURF_MESH_R/ + & MSHBLK(3, 4*NVMAX), ! block to store all surface meshes + & RV1MSH(3,NVMAX), ! mesh h.v. vortex left points + & RV2MSH(3,NVMAX), ! mesh h.v. vortex right points + & RVMSH(3,NVMAX), ! mesh h.v. vortex center points + & RCMSH(3,NVMAX) ! mesh h.v. control points + + LOGICAL LSURFMSH + LOGICAL LMESHFLAT + COMMON /SURF_MESH_L/ + & LSURFMSH(NFMAX), ! T if surface uses a mesh cordinates for its geometry + & LMESHFLAT(NFMAX) ! T if the surface's mesh should be flattened when computing vorticies and control points + + +C !!--- end added variables for python geometry manipulation --- COMMON /STRP_I/ & NSURFS(NSMAX), ! index of surface which contains this strip @@ -518,6 +540,7 @@ C !!--- end added variables for python geometry minipulation --- REAL(kind=avl_real) RLE2, CHORD2 REAL(kind=avl_real) WSTRIP REAL(kind=avl_real) TANLE, TANTE + REAL(kind=avl_real) GINCSTRIP REAL(kind=avl_real) CLCD REAL(kind=avl_real) SAXFR REAL(kind=avl_real) ESS @@ -553,6 +576,7 @@ C !!--- end added variables for python geometry minipulation --- & RLE2(3,NSMAX), CHORD2(NSMAX), ! strip right end LE point, chord & WSTRIP(NSMAX), ! strip y-z width & TANLE(NSMAX), TANTE(NSMAX), ! strip LE,TE sweep slopes + & GINCSTRIP(NSMAX), ! strip geometric incidence angle & CLCD(NUMAX,NSMAX), ! strip viscous polar & SAXFR, ! x/c of spanwise axis for Vperp def & ESS(3,NSMAX), ! spanwise unit vector for Vperp def @@ -610,6 +634,7 @@ C REAL(kind=avl_real) RS REAL(kind=avl_real) RL,RADL REAL(kind=avl_real) DXV + REAL(kind=avl_real) DXSTRPV REAL(kind=avl_real) CHORDV REAL(kind=avl_real) SLOPEV REAL(kind=avl_real) SLOPEC @@ -654,6 +679,7 @@ C & RS(3,NVMAX), ! h.v. source points & RL(3,NLMAX),RADL(NLMAX), ! source line node points, body radius & DXV(NVMAX), ! chord of element + & DXSTRPV(NVMAX), ! chord of element projected onto chord of element-containing strip & CHORDV(NVMAX), ! chord of element-containing strip & SLOPEV(NVMAX), ! camber slopes at h.v. bound leg & SLOPEC(NVMAX), ! camber slopes at c.p. diff --git a/src/includes/AVL_ad_seeds.inc b/src/includes/AVL_ad_seeds.inc index 1900531..312a556 100644 --- a/src/includes/AVL_ad_seeds.inc +++ b/src/includes/AVL_ad_seeds.inc @@ -354,6 +354,8 @@ C real(kind=avl_real) XBOD_DIFF real(kind=avl_real) YBOD_DIFF real(kind=avl_real) TBOD_DIFF + real(kind=avl_real) XBOD_R_DIFF + real(kind=avl_real) YBOD_R_DIFF COMMON /BODY_GEOM_R_DIFF/ & XYZSCAL_B_DIFF(3, NBMAX), & XYZTRAN_B_DIFF(3, NBMAX), @@ -362,7 +364,9 @@ C & BSPACE_DIFF(NBMAX), & XBOD_DIFF(IBX, NBMAX), & YBOD_DIFF(IBX, NBMAX), - & TBOD_DIFF(IBX, NBMAX) + & TBOD_DIFF(IBX, NBMAX), + & XBOD_R_DIFF(IBX, NBMAX), + & YBOD_R_DIFF(IBX, NBMAX) C real(kind=avl_real) XYZSCAL_DIFF real(kind=avl_real) XYZTRAN_DIFF @@ -374,8 +378,8 @@ C real(kind=avl_real) XYZLES_DIFF real(kind=avl_real) XSEC_DIFF real(kind=avl_real) YSEC_DIFF - real(kind=avl_real) XFMIN_DIFF - real(kind=avl_real) XFMAX_DIFF + real(kind=avl_real) XFMIN_R_DIFF + real(kind=avl_real) XFMAX_R_DIFF real(kind=avl_real) XLASEC_DIFF real(kind=avl_real) ZLASEC_DIFF real(kind=avl_real) XUASEC_DIFF @@ -405,8 +409,8 @@ C & XYZLES_DIFF(3, NSMAX, NFMAX), & XSEC_DIFF(IBX, NSMAX, NFMAX), & YSEC_DIFF(IBX, NSMAX, NFMAX), - & XFMIN_DIFF(NSMAX, NFMAX), - & XFMAX_DIFF(NSMAX, NFMAX), + & XFMIN_R_DIFF(NSMAX, NFMAX), + & XFMAX_R_DIFF(NSMAX, NFMAX), & XLASEC_DIFF(IBX, NSMAX, NFMAX), & ZLASEC_DIFF(IBX, NSMAX, NFMAX), & XUASEC_DIFF(IBX, NSMAX, NFMAX), @@ -425,6 +429,10 @@ C & GAIND_DIFF(ICONX, NSMAX, NFMAX), & REFLD_DIFF(ICONX, NSMAX, NFMAX), & GAING_DIFF(ICONX, NSMAX, NFMax) +C + real(kind=avl_real) MSHBLK_DIFF + COMMON /SURF_MESH_R_DIFF/ + & MSHBLK_DIFF(3, 4*NVMAX) C real(kind=avl_real) RLE_DIFF real(kind=avl_real) CHORD_DIFF @@ -435,6 +443,7 @@ C real(kind=avl_real) WSTRIP_DIFF real(kind=avl_real) TANLE_DIFF real(kind=avl_real) TANTE_DIFF + real(kind=avl_real) GINCSTRIP_DIFF real(kind=avl_real) CLCD_DIFF real(kind=avl_real) SAXFR_DIFF real(kind=avl_real) ESS_DIFF @@ -505,6 +514,7 @@ C & WSTRIP_DIFF(NSMAX), & TANLE_DIFF(NSMAX), & TANTE_DIFF(NSMAX), + & GINCSTRIP_DIFF(NSMAX), & CLCD_DIFF(NUMAX,NSMAX), & SAXFR_DIFF, & ESS_DIFF(3,NSMAX), @@ -574,6 +584,7 @@ C real(kind=avl_real) RL_DIFF real(kind=avl_real) RADL_DIFF real(kind=avl_real) DXV_DIFF + real(kind=avl_real) DXSTRPV_DIFF real(kind=avl_real) CHORDV_DIFF real(kind=avl_real) SLOPEV_DIFF real(kind=avl_real) SLOPEC_DIFF @@ -619,6 +630,7 @@ C & RL_DIFF(3,NLMAX), & RADL_DIFF(NLMAX), & DXV_DIFF(NVMAX), + & DXSTRPV_DIFF(NVMAX), & CHORDV_DIFF(NVMAX), & SLOPEV_DIFF(NVMAX), & SLOPEC_DIFF(NVMAX), diff --git a/src/sgutil.f b/src/sgutil.f index b869f5e..8ab763f 100644 --- a/src/sgutil.f +++ b/src/sgutil.f @@ -296,6 +296,24 @@ SUBROUTINE SPACER (N,PSPACE,X) SUBROUTINE CSPACER(NVC,CSPACE,CLAF, XPT,XVR,XSR,XCP) + ! This is a extremely important funciton that is not + ! documented for some reason. + ! Inputs: + ! NVC: NUMBER OF DESIRED POINTS IN ARRAY + ! CSPACE: SPACING PARAMETER (-3<=PSPACE<=3). + ! DEFINES POINT DISTRIBUTION + ! TO BE USED AS FOLLOWS: + ! PSPACE = 0 : EQUAL SPACING + ! PSPACE = 1 : COSINE SPACING. + ! PSPACE = 2 : SINE SPACING + ! (CONCENTRATING POINTS NEAR 0). + ! PSPACE = 3 : EQUAL SPACING. + ! CLAF: CL alfa (needed to determine control point location) + ! Outputs: + ! XPT: Array of panel leading edge x-locations + ! XVR: Array of vortex x-locations + ! XSR: Array of source x-locations + ! XCP: Array of control point x-locations REAL XPT(*), XVR(*), XSR(*), XCP(*) C PI = 4.0*ATAN(1.0) @@ -318,17 +336,21 @@ SUBROUTINE CSPACER(NVC,CSPACE,CLAF, XPT,XVR,XSR,XCP) ENDIF C C---- cosine chordwise spacing + ! Each of these provides a quarter panel chord offset for cosine, + ! sine, and uniform spacing respectively. DTH1 = PI/FLOAT(4*NVC + 2) DTH2 = 0.5*PI/FLOAT(4*NVC + 1) DXC0 = 1.0/FLOAT(4*NVC) C DO IVC = 1, NVC C------ uniform - XC0 = INT(4*IVC - 4) * DXC0 + XC0 = INT(4*IVC - 4) * DXC0 ! eqv (IVC-1)/NVC XPT0 = XC0 - XVR0 = XC0 + DXC0 - XSR0 = XC0 + 2.0*DXC0 - XCP0 = XC0 + DXC0 + 2.0*DXC0*CLAF + XVR0 = XC0 + DXC0 ! quarter-chord + XSR0 = XC0 + 2.0*DXC0 ! half-chord + XCP0 = XC0 + DXC0 + 2.0*DXC0*CLAF ! quarter-chord + half-chord*claf + ! Note: claf is a scaling factor so typically claf = 1 and the control point + ! is at the three-quarter chord position of the panel C C------ cosine TH1 = INT(4*IVC - 3) * DTH1