diff --git a/.github/workflows/pylint.yml b/.github/workflows/pylint.yml new file mode 100644 index 0000000..c73e032 --- /dev/null +++ b/.github/workflows/pylint.yml @@ -0,0 +1,23 @@ +name: Pylint + +on: [push] + +jobs: + build: + runs-on: ubuntu-latest + strategy: + matrix: + python-version: ["3.8", "3.9", "3.10"] + steps: + - uses: actions/checkout@v4 + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v3 + with: + python-version: ${{ matrix.python-version }} + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install pylint + - name: Analysing the code with pylint + run: | + pylint $(git ls-files '*.py') diff --git a/.gitignore b/.gitignore old mode 100644 new mode 100755 diff --git a/GeoMACH/BSE/BSEmodel.py b/GeoMACH/BSE/BSEmodel.py old mode 100644 new mode 100755 index 0f9974d..4bb00fb --- a/GeoMACH/BSE/BSEmodel.py +++ b/GeoMACH/BSE/BSEmodel.py @@ -459,10 +459,12 @@ def compute_projection(self, name, pts, surf_pts=None, ndim=1): npts = pts.shape[0] nsurf_pts = surf_pts.shape[0] + nrefine = 10 surfs, ind_u, ind_v \ = BSElib.computeproj(npts, nsurf_pts, size['cp_str'], size['pt_str'], - num['surf'], num['group'], surf_pts, + num['surf'], num['group'], nrefine, + surf_pts, str_indices['cp'], str_indices['pt'], topo['surf_group'], bspline['order'], bspline['num_cp'], bspline['num_pt'], diff --git a/GeoMACH/BSE/BSEvec.py b/GeoMACH/BSE/BSEvec.py old mode 100644 new mode 100755 diff --git a/GeoMACH/BSE/__init__.py b/GeoMACH/BSE/__init__.py old mode 100644 new mode 100755 diff --git a/GeoMACH/BSE/test.py b/GeoMACH/BSE/test.py old mode 100644 new mode 100755 index 4bf4981..78ce35b --- a/GeoMACH/BSE/test.py +++ b/GeoMACH/BSE/test.py @@ -55,6 +55,7 @@ def cube(nx, ny, nz, rx, ry, rz): surf.compute_projection('test', pts0, ndim=3) surf.apply_jacobian('test', 'd(test)/d(cp_str)', 'cp_str') pts = numpy.array(surf.vec['test'].array) +print pts #surf.vec['df'].export_tec_scatter() surf.vec['pt_str'].export_tec_str() #surf.vec['pt_str'].export_STL() diff --git a/GeoMACH/PGM/__init__.py b/GeoMACH/PGM/__init__.py old mode 100644 new mode 100755 diff --git a/GeoMACH/PGM/airfoils/n64206.dat b/GeoMACH/PGM/airfoils/n64206.dat old mode 100644 new mode 100755 diff --git a/GeoMACH/PGM/airfoils/rae2822.dat b/GeoMACH/PGM/airfoils/rae2822.dat old mode 100644 new mode 100755 diff --git a/GeoMACH/PGM/components/PGMbody.py b/GeoMACH/PGM/components/PGMbody.py old mode 100644 new mode 100755 diff --git a/GeoMACH/PGM/components/PGMcone.py b/GeoMACH/PGM/components/PGMcone.py old mode 100644 new mode 100755 diff --git a/GeoMACH/PGM/components/PGMinterpolant.py b/GeoMACH/PGM/components/PGMinterpolant.py old mode 100644 new mode 100755 index 2dc32da..b8c1e69 --- a/GeoMACH/PGM/components/PGMinterpolant.py +++ b/GeoMACH/PGM/components/PGMinterpolant.py @@ -95,7 +95,7 @@ def compute_normals(self): self._bse.apply_jacobian(qname, 'd(' + qname + ')/d(cp_str)', 'cp_str') cross = numpy.cross(self._bse.vec[name + '_du'].array, self._bse.vec[name + '_dv'].array) - norms = numpy.linalg.norm(cross, axis=1) + norms = numpy.sqrt(numpy.sum(cross**2,axis=1)) for ind in xrange(3): cross[:, ind] /= norms diff --git a/GeoMACH/PGM/components/PGMjunction.py b/GeoMACH/PGM/components/PGMjunction.py old mode 100644 new mode 100755 diff --git a/GeoMACH/PGM/components/PGMprimitive.py b/GeoMACH/PGM/components/PGMprimitive.py old mode 100644 new mode 100755 diff --git a/GeoMACH/PGM/components/PGMshell.py b/GeoMACH/PGM/components/PGMshell.py old mode 100644 new mode 100755 diff --git a/GeoMACH/PGM/components/PGMtip.py b/GeoMACH/PGM/components/PGMtip.py old mode 100644 new mode 100755 diff --git a/GeoMACH/PGM/components/PGMwing.py b/GeoMACH/PGM/components/PGMwing.py old mode 100644 new mode 100755 diff --git a/GeoMACH/PGM/components/__init__.py b/GeoMACH/PGM/components/__init__.py old mode 100644 new mode 100755 diff --git a/GeoMACH/PGM/core/MACHconfiguration.py b/GeoMACH/PGM/core/MACHconfiguration.py old mode 100644 new mode 100755 index d4e7cfd..f2f9fef --- a/GeoMACH/PGM/core/MACHconfiguration.py +++ b/GeoMACH/PGM/core/MACHconfiguration.py @@ -15,6 +15,7 @@ def __init__(self): self.diff = OrderedDict() self.jacobians = OrderedDict() self.updated = {} + self.debug = False def addPointSet(self, points, pt_name, origConfig=True, **kwargs): bse = self._bse @@ -26,6 +27,13 @@ def addPointSet(self, points, pt_name, origConfig=True, **kwargs): self.jacobians[pt_name] = bse.jac['d(' + pt_name + ')/d(cp_str)'] self.diff[pt_name] = points - self.jacobians[pt_name].dot(bse.vec['cp_str'].array) + ### FOR DEBUGGING (Print projected points) + if self.debug is True: + bse.apply_jacobian(pt_name, 'd(' + pt_name + ')/d(cp_str)', 'cp_str') + bse.vec[pt_name].export_tec_scatter('projected_points.dat') + bse.vec[pt_name].array[:, :] = points + bse.vec[pt_name].export_tec_scatter('CFD_surf_points.dat') + self.updated[pt_name] = False def setDesignVars(self, dv_dict): diff --git a/GeoMACH/PGM/core/PGMcomponent.py b/GeoMACH/PGM/core/PGMcomponent.py old mode 100644 new mode 100755 index c1cee9d..3e3fe2d --- a/GeoMACH/PGM/core/PGMcomponent.py +++ b/GeoMACH/PGM/core/PGMcomponent.py @@ -42,6 +42,8 @@ def __init__(self): self._num_surf = {} self._shapes = {} + self.funcs = {} + def initialize_props(self): """ Adds the *X*, *Y*, and *Z* shape properties """ props = self.props diff --git a/GeoMACH/PGM/core/PGMconfiguration.py b/GeoMACH/PGM/core/PGMconfiguration.py old mode 100644 new mode 100755 diff --git a/GeoMACH/PGM/core/PGMdv.py b/GeoMACH/PGM/core/PGMdv.py old mode 100644 new mode 100755 diff --git a/GeoMACH/PGM/core/PGMface.py b/GeoMACH/PGM/core/PGMface.py old mode 100644 new mode 100755 diff --git a/GeoMACH/PGM/core/PGMobject.py b/GeoMACH/PGM/core/PGMobject.py old mode 100644 new mode 100755 diff --git a/GeoMACH/PGM/core/PGMparameter.py b/GeoMACH/PGM/core/PGMparameter.py old mode 100644 new mode 100755 diff --git a/GeoMACH/PGM/core/PGMproperty.py b/GeoMACH/PGM/core/PGMproperty.py old mode 100644 new mode 100755 diff --git a/GeoMACH/PGM/core/PGMsurf.py b/GeoMACH/PGM/core/PGMsurf.py old mode 100644 new mode 100755 diff --git a/GeoMACH/PGM/core/PGMvec.py b/GeoMACH/PGM/core/PGMvec.py old mode 100644 new mode 100755 diff --git a/GeoMACH/PGM/core/__init__.py b/GeoMACH/PGM/core/__init__.py old mode 100644 new mode 100755 diff --git a/GeoMACH/PSM/BDFwriter.py b/GeoMACH/PSM/BDFwriter.py old mode 100644 new mode 100755 diff --git a/GeoMACH/PSM/QUAD.py b/GeoMACH/PSM/QUAD.py old mode 100644 new mode 100755 diff --git a/GeoMACH/PSM/__init__.py b/GeoMACH/PSM/__init__.py old mode 100644 new mode 100755 diff --git a/GeoMACH/PSM/airframe.py b/GeoMACH/PSM/airframe.py old mode 100644 new mode 100755 diff --git a/GeoMACH/__init__.py b/GeoMACH/__init__.py old mode 100644 new mode 100755 diff --git a/MANIFEST.in b/MANIFEST.in old mode 100644 new mode 100755 diff --git a/README.txt b/README.txt old mode 100644 new mode 100755 index 126f8c9..99ae1c0 --- a/README.txt +++ b/README.txt @@ -28,3 +28,21 @@ Installing and getting started (linux or Mac) 5. Run the example: $ cd [GeoMACH-top]/examples $ python conventional.py + +Installing on a local directory (linux or Mac, for clusters) +--------------------------------------------- +[GeoMACH-top] = path to the top-level GeoMACH directory + +1. Install (or load) numpy, scipy + +2. Install GeoMACH: + $ cd [GeoMACH-top] + $ python setup.py build --fcompiler=intelem (you may use another compiler) + $ python setup.py install --user + +3. You have to manually copy all airfoils to the installation folder + $ cp -r GeoMACH/PGM/airfoils ~/.local/lib/python2.7/site-packages/GeoMACH-0.1-py2.7-linux-x86_64.egg/GeoMACH/PGM/ + +4. Run the example: + $ cd examples + $ python conventional.py diff --git a/docs/Makefile b/docs/Makefile old mode 100644 new mode 100755 diff --git a/docs/conf.py b/docs/conf.py old mode 100644 new mode 100755 diff --git a/docs/index.rst b/docs/index.rst old mode 100644 new mode 100755 diff --git a/docs/src/Conventional.rst b/docs/src/Conventional.rst old mode 100644 new mode 100755 diff --git a/docs/src/GeoMACH.BSE.rst b/docs/src/GeoMACH.BSE.rst old mode 100644 new mode 100755 diff --git a/docs/src/GeoMACH.PGM.components.rst b/docs/src/GeoMACH.PGM.components.rst old mode 100644 new mode 100755 diff --git a/docs/src/GeoMACH.PGM.core.rst b/docs/src/GeoMACH.PGM.core.rst old mode 100644 new mode 100755 diff --git a/docs/src/GeoMACH.PGM.rst b/docs/src/GeoMACH.PGM.rst old mode 100644 new mode 100755 diff --git a/docs/src/GeoMACH.rst b/docs/src/GeoMACH.rst old mode 100644 new mode 100755 diff --git a/docs/src/Images/UML2.svg b/docs/src/Images/UML2.svg old mode 100644 new mode 100755 diff --git a/docs/src/Images/body.png b/docs/src/Images/body.png old mode 100644 new mode 100755 diff --git a/docs/src/Images/bsplines.png b/docs/src/Images/bsplines.png old mode 100644 new mode 100755 diff --git a/docs/src/Images/cone.png b/docs/src/Images/cone.png old mode 100644 new mode 100755 diff --git a/docs/src/Images/conventional.jpg b/docs/src/Images/conventional.jpg old mode 100644 new mode 100755 diff --git a/docs/src/Images/cp.png b/docs/src/Images/cp.png old mode 100644 new mode 100755 diff --git a/docs/src/Images/d8.jpg b/docs/src/Images/d8.jpg old mode 100644 new mode 100755 diff --git a/docs/src/Images/df.png b/docs/src/Images/df.png old mode 100644 new mode 100755 diff --git a/docs/src/Images/geom1.png b/docs/src/Images/geom1.png old mode 100644 new mode 100755 diff --git a/docs/src/Images/hwb.png b/docs/src/Images/hwb.png old mode 100644 new mode 100755 diff --git a/docs/src/Images/joinedwing.jpg b/docs/src/Images/joinedwing.jpg old mode 100644 new mode 100755 diff --git a/docs/src/Images/junction.jpg b/docs/src/Images/junction.jpg old mode 100644 new mode 100755 diff --git a/docs/src/Images/junction.png b/docs/src/Images/junction.png old mode 100644 new mode 100755 diff --git a/docs/src/Images/junction1.jpg b/docs/src/Images/junction1.jpg old mode 100644 new mode 100755 diff --git a/docs/src/Images/primitives.png b/docs/src/Images/primitives.png old mode 100644 new mode 100755 diff --git a/docs/src/Images/primitives0.png b/docs/src/Images/primitives0.png old mode 100644 new mode 100755 diff --git a/docs/src/Images/pt.png b/docs/src/Images/pt.png old mode 100644 new mode 100755 diff --git a/docs/src/Images/shell.png b/docs/src/Images/shell.png old mode 100644 new mode 100755 diff --git a/docs/src/Images/strutbraced.jpg b/docs/src/Images/strutbraced.jpg old mode 100644 new mode 100755 diff --git a/docs/src/Images/supersonic.png b/docs/src/Images/supersonic.png old mode 100644 new mode 100755 diff --git a/docs/src/Images/tip.png b/docs/src/Images/tip.png old mode 100644 new mode 100755 diff --git a/docs/src/Images/wing.png b/docs/src/Images/wing.png old mode 100644 new mode 100755 diff --git a/docs/src/Tutorials.rst b/docs/src/Tutorials.rst old mode 100644 new mode 100755 diff --git a/examples/RC_aircraft_tutorial.py b/examples/RC_aircraft_tutorial.py old mode 100644 new mode 100755 diff --git a/examples/conventional.py b/examples/conventional.py old mode 100644 new mode 100755 diff --git a/examples/supersonic.py b/examples/supersonic.py old mode 100644 new mode 100755 diff --git a/examples/trussbraced_full.py b/examples/trussbraced_full.py old mode 100644 new mode 100755 diff --git a/examples/trussbraced_wingstrut.py b/examples/trussbraced_wingstrut.py old mode 100644 new mode 100755 diff --git a/examples/wing.py b/examples/wing.py old mode 100644 new mode 100755 index b618bcd..6a0c731 --- a/examples/wing.py +++ b/examples/wing.py @@ -62,6 +62,7 @@ def meshStructure(self, res, filename): bse = pgm.initialize() pgm.comps['wing'].set_airfoil('rae2822.dat') bse.vec['pt_str'].export_tec_str() - + bse.vec['cp_str'].export_IGES() + exit() for num in [20]:#[5, 10, 20, 50, 100]: pgm.meshStructure(num, 'wing_'+str(num)) diff --git a/setup.cfg b/setup.cfg old mode 100644 new mode 100755 diff --git a/setup.py b/setup.py old mode 100644 new mode 100755 diff --git a/src/BSE/bspline_basis.f90 b/src/BSE/bspline_basis.f90 old mode 100644 new mode 100755 diff --git a/src/BSE/bspline_knot.f90 b/src/BSE/bspline_knot.f90 old mode 100644 new mode 100755 diff --git a/src/BSE/bspline_param.f90 b/src/BSE/bspline_param.f90 old mode 100644 new mode 100755 diff --git a/src/BSE/compute_bs_jacobian.f90 b/src/BSE/compute_bs_jacobian.f90 old mode 100644 new mode 100755 diff --git a/src/BSE/compute_cp_jacobian.f90 b/src/BSE/compute_cp_jacobian.f90 old mode 100644 new mode 100755 diff --git a/src/BSE/compute_df_jacobian.f90 b/src/BSE/compute_df_jacobian.f90 old mode 100644 new mode 100755 diff --git a/src/BSE/compute_in_jacobian.f90 b/src/BSE/compute_in_jacobian.f90 old mode 100644 new mode 100755 diff --git a/src/BSE/compute_indices.f90 b/src/BSE/compute_indices.f90 old mode 100644 new mode 100755 diff --git a/src/BSE/compute_projection.f90 b/src/BSE/compute_projection.f90 old mode 100644 new mode 100755 index 6d53cc8..465a7ad --- a/src/BSE/compute_projection.f90 +++ b/src/BSE/compute_projection.f90 @@ -1,5 +1,5 @@ subroutine computeproj(npt_proj, nsurf_proj, & - ncp_str, npt_str, nsurf, ngroup, & + ncp_str, npt_str, nsurf, ngroup, nrefine, & surf_proj, str_indices_cp, str_indices_pt, & surf_group, order, num_cp, num_pt, & cp_str, pt_str, pt_proj, & @@ -8,7 +8,7 @@ subroutine computeproj(npt_proj, nsurf_proj, & implicit none !Fortran-python interface directives - !f2py intent(in) npt_proj, nsurf_proj, ncp_str, npt_str, nsurf, ngroup, surf_proj, str_indices_cp, str_indices_pt, surf_group, order, num_cp, num_pt, cp_str, pt_str, pt_proj + !f2py intent(in) npt_proj, nsurf_proj, ncp_str, npt_str, nsurf, ngroup, nrefine, surf_proj, str_indices_cp, str_indices_pt, surf_group, order, num_cp, num_pt, cp_str, pt_str, pt_proj !f2py intent(out) min_s, min_tu, min_tv !f2py depend(nsurf_proj) surf_proj !f2py depend(nsurf) str_indices_cp, str_indices_pt, surf_group @@ -20,7 +20,7 @@ subroutine computeproj(npt_proj, nsurf_proj, & !Input integer, intent(in) :: npt_proj, nsurf_proj - integer, intent(in) :: ncp_str, npt_str, nsurf, ngroup + integer, intent(in) :: ncp_str, npt_str, nsurf, ngroup, nrefine integer, intent(in) :: surf_proj(nsurf_proj) integer, intent(in) :: str_indices_cp(nsurf, 2) integer, intent(in) :: str_indices_pt(nsurf, 2) @@ -40,16 +40,20 @@ subroutine computeproj(npt_proj, nsurf_proj, & integer isurf_proj, isurf, ipt_proj integer offset_cp, offset_pt, offset_u, offset_v integer ku, kv, mu, mv, nu, nv - integer u, v, min_u, min_v - integer lu, lv, hu, hv, k, counter - double precision pt0(3), pt(3), min_d_surf, d - double precision P(3), Pu(3), Pv(3), Puu(3), Puv(3), Pvv(3) - double precision f(3), x(2), g(2), H(2, 2), det - double precision norm_g, norm_dx, W(2, 2), dx(2) + integer u, v, uu, vv + integer u1, u2, v1, v2 + integer lu, lv, hu, hv + integer min_s_local + double precision tu, tv + double precision min_d_local, min_tu_local, min_tv_local + double precision pt0(3), pt(3), d + double precision P(3), P00(3), P01(3), P10(3), P11(3) + double precision x(2), x0(2) + logical fail + double precision normal0(3), normal1(3), normal2(3), normal3(3), normal4(3) double precision, allocatable, dimension(:) :: knot_u, knot_v double precision, allocatable, dimension(:) :: param_u, param_v - double precision, allocatable, dimension(:) :: bu0, bu1, bu2 - double precision, allocatable, dimension(:) :: bv0, bv1, bv2 + double precision, allocatable, dimension(:) :: bu0, bv0 double precision, allocatable, dimension(:, :, :) :: cp min_d(:) = 1e10 @@ -82,11 +86,7 @@ subroutine computeproj(npt_proj, nsurf_proj, & call paramuni(kv+mv, mv, nv, knot_v, param_v) allocate(bu0(ku)) - allocate(bu1(ku)) - allocate(bu2(ku)) allocate(bv0(kv)) - allocate(bv1(kv)) - allocate(bv2(kv)) allocate(cp(mu, mv, 3)) do u = 1, mu @@ -99,112 +99,130 @@ subroutine computeproj(npt_proj, nsurf_proj, & do ipt_proj = 1, npt_proj pt0 = pt_proj(ipt_proj, :) - min_d_surf = 1e10 - min_u = 1 - min_v = 1 + min_d_local = 1e10 + min_s_local = 1 + min_tu_local = 0.0 + min_tv_local = 0.0 do u = 1, nu ! ceiling(nu/100.0) do v = 1, nv ! ceiling(nv/100.0) pt = pt_str(offset_pt + (v-1)*nu + u, :) d = abs(dot_product(pt - pt0, pt - pt0)) - if (d .lt. min_d_surf) then - min_d_surf = d - min_u = u - min_v = v + if (d .lt. min_d_local) then + min_d_local = d + min_s_local = isurf + min_tu_local = param_u(u) + min_tv_local = param_v(v) end if end do end do - if (min_d_surf .lt. min_d(ipt_proj)) then - min_d(ipt_proj) = min_d_surf - min_s(ipt_proj) = isurf - min_tu(ipt_proj) = param_u(min_u) - min_tv(ipt_proj) = param_v(min_v) + if (min_d_local .lt. min_d(ipt_proj)) then + min_d(ipt_proj) = min_d_local + min_s(ipt_proj) = min_s_local + min_tu(ipt_proj) = min_tu_local + min_tv(ipt_proj) = min_tv_local end if - x(1) = param_u(min_u) - x(2) = param_v(min_v) - do counter = 0, 40 - call basis0(ku, ku+mu, x(1), knot_u, bu0, offset_u) - call basis1(ku, ku+mu, x(1), knot_u, bu1, offset_u) - call basis2(ku, ku+mu, x(1), knot_u, bu2, offset_u) - call basis0(kv, kv+mv, x(2), knot_v, bv0, offset_v) - call basis1(kv, kv+mv, x(2), knot_v, bv1, offset_v) - call basis2(kv, kv+mv, x(2), knot_v, bv2, offset_v) - P(:) = 0.0 - Pu(:) = 0.0 - Pv(:) = 0.0 - Puu(:) = 0.0 - Puv(:) = 0.0 - Pvv(:) = 0.0 - do lu = 1, ku - hu = lu + offset_u - do lv = 1, kv - hv = lv + offset_v - P = P + bu0(lu) * bv0(lv) * cp(hu, hv, :) - Pu = Pu + bu1(lu) * bv0(lv) * cp(hu, hv, :) - Pv = Pv + bu0(lu) * bv1(lv) * cp(hu, hv, :) - Puu = Puu + bu2(lu) * bv0(lv) * cp(hu, hv, :) - Puv = Puv + bu1(lu) * bv1(lv) * cp(hu, hv, :) - Pvv = Pvv + bu0(lu) * bv2(lv) * cp(hu, hv, :) - end do - end do + x0(1) = min_tu_local + x0(2) = min_tv_local + call newtonprojection(ku, kv, mu, mv, x0, pt0, cp, knot_u, knot_v, x, pt, fail) - f = P - pt0 - g(1) = 2 * dot_product(f, Pu) - g(2) = 2 * dot_product(f, Pv) - H(1, 1) = 2 * dot_product(Pu, Pu) + 2 * dot_product(f, Puu) - H(1, 2) = 2 * dot_product(Pu, Pv) + 2 * dot_product(f, Puv) - H(2, 2) = 2 * dot_product(Pv, Pv) + 2 * dot_product(f, Pvv) - H(2, 1) = H(1, 2) - do k = 1, 2 - if (((x(k) .eq. 0) .and. (g(k) .gt. 0)) .or. & - ((x(k) .eq. 1) .and. (g(k) .lt. 0))) then - g(k) = 0.0 - H(1, 2) = 0.0 - H(2, 1) = 0.0 - H(k, k) = 1.0 - end if - end do - det = H(1, 1) * H(2, 2) - H(1, 2) * H(2, 1) - W(1, 1) = H(2, 2) / det - W(1, 2) = -H(1, 2) / det - W(2, 2) = H(1, 1) / det - W(2, 1) = W(1, 2) - norm_g = sqrt(dot_product(g, g)) - dx(1) = -dot_product(W(1, :), g) - dx(2) = -dot_product(W(2, :), g) - do k = 1, 2 - if (x(k) + dx(k) .lt. 0) then - dx(k) = -x(k) - else if (x(k) + dx(k) .gt. 1) then - dx(k) = 1 - x(k) - end if - end do - norm_dx = sqrt(dot_product(dx, dx)) - - ! print *, counter, norm - if ((norm_g .lt. 1e-13) .or. (norm_dx .lt. 1e-13)) then - exit - end if - x = x + dx - end do - - call basis0(ku, ku+mu, x(1), knot_u, bu0, offset_u) - call basis0(kv, kv+mv, x(2), knot_v, bv0, offset_v) - P(:) = 0.0 - do lu = 1, ku - hu = lu + offset_u - do lv = 1, kv - hv = lv + offset_v - P = P + bu0(lu) * bv0(lv) * cp(hu, hv, :) - end do - end do - d = abs(dot_product(P - pt0, P - pt0)) + d = abs(dot_product(pt - pt0, pt - pt0)) if (d .lt. min_d(ipt_proj)) then min_d(ipt_proj) = d min_s(ipt_proj) = isurf min_tu(ipt_proj) = x(1) min_tv(ipt_proj) = x(2) end if + + !Check if Newton Search converged + if (fail) then !We will refine a cell to get a better starting point for the Newton Search + do u = 1, nu-1 !Looping over each surface point + u1 = u + u2 = u + 1 + do v = 1, nv-1 + v1 = v + v2 = v + 1 + + !Get the corner points of the cell + P00 = pt_str(offset_pt + (v1-1)*nu + u1, :) + P10 = pt_str(offset_pt + (v1-1)*nu + u2, :) + P01 = pt_str(offset_pt + (v2-1)*nu + u1, :) + P11 = pt_str(offset_pt + (v2-1)*nu + u2, :) + call cross_product(P11 - P00, P01 - P10, normal0) !Find a vector normal to the cell + + !Cross product between cell edges and the CFD point position with respect to point P00 + call cross_product(pt0 - P00, P01 - P00, normal1) !u=0 + call cross_product(pt0 - P01, P11 - P01, normal2) !v=1 + call cross_product(pt0 - P11, P10 - P11, normal3) !u=1 + call cross_product(pt0 - P10, P00 - P10, normal4) !v=0 + + !Identify if the point is above the cell + if (& + (dot_product(normal1, normal0) .ge. 0) .and. & + (dot_product(normal2, normal0) .ge. 0) .and. & + (dot_product(normal3, normal0) .ge. 0) .and. & + (dot_product(normal4, normal0) .ge. 0)) then + min_d_local = 1e10 + min_s_local = 1 + min_tu_local = 0.0 + min_tv_local = 0.0 + !Loop over refined points + do uu = 1, nrefine + tu = param_u(u1) + (param_u(u2) - param_u(u1)) * (uu-1) / (nrefine-1) + do vv = 1, nrefine + tv = param_v(v1) + (param_v(v2) - param_v(v1)) * (vv-1) / (nrefine-1) + + !Compute physical coordinates of refined point + call basis0(ku, ku+mu, tu, knot_u, bu0, offset_u) + call basis0(kv, kv+mv, tv, knot_v, bv0, offset_v) + P(:) = 0.0 + do lu = 1, ku + hu = lu + offset_u + do lv = 1, kv + hv = lv + offset_v + P = P + bu0(lu) * bv0(lv) * cp(hu, hv, :) + end do + end do + pt(:) = P(:) + + !Compare refined point to the CFD mesh point + d = abs(dot_product(pt - pt0, pt - pt0)) + if (d .lt. min_d_local) then + min_d_local = d + min_s_local = isurf + min_tu_local = tu + min_tv_local = tv + end if + end do + end do + + !Compare best candidate with best candidate prior to the refinement + if (min_d_local .lt. min_d(ipt_proj)) then + min_d(ipt_proj) = min_d_local + min_s(ipt_proj) = min_s_local + min_tu(ipt_proj) = min_tu_local + min_tv(ipt_proj) = min_tv_local + end if + + !Call newton search to find u and v values + x0(1) = min_tu_local + x0(2) = min_tv_local + call newtonprojection(ku, kv, mu, mv, x0, pt0, cp, knot_u, knot_v, x, pt, fail) + + d = abs(dot_product(pt - pt0, pt - pt0)) + if (d .lt. min_d(ipt_proj)) then + min_d(ipt_proj) = d + min_s(ipt_proj) = isurf + min_tu(ipt_proj) = x(1) + min_tv(ipt_proj) = x(2) + end if + + end if + + end do + end do + end if + end do deallocate(knot_u) @@ -214,13 +232,128 @@ subroutine computeproj(npt_proj, nsurf_proj, & deallocate(param_v) deallocate(bu0) - deallocate(bu1) - deallocate(bu2) deallocate(bv0) - deallocate(bv1) - deallocate(bv2) deallocate(cp) end do end subroutine computeproj + + + +subroutine newtonprojection(ku, kv, mu, mv, & + x0, pt0, cp, knot_u, knot_v, x, pt, fail) + + implicit none + + !Input + integer, intent(in) :: ku, kv, mu, mv + double precision, intent(in) :: x0(2), pt0(3) + double precision, intent(in) :: cp(mu, mv, 3) + double precision, intent(in) :: knot_u(ku + mu), knot_v(kv + mv) + + !Output + double precision, intent(out) :: x(2), pt(3) + logical, intent(out) :: fail + + !Working + integer offset_u, offset_v + double precision bu0(ku), bu1(ku), bu2(ku) + double precision bv0(kv), bv1(kv), bv2(kv) + double precision P(3), Pu(3), Pv(3), Puu(3), Puv(3), Pvv(3) + double precision f(3), g(2), H(2, 2), det + double precision norm_g, norm_dx, norm_f, W(2, 2), dx(2) + integer lu, lv, hu, hv, k, counter + + x(:) = x0(:) + fail = .True. + + do counter = 0, 40 + call basis0(ku, ku+mu, x(1), knot_u, bu0, offset_u) + call basis1(ku, ku+mu, x(1), knot_u, bu1, offset_u) + call basis2(ku, ku+mu, x(1), knot_u, bu2, offset_u) + call basis0(kv, kv+mv, x(2), knot_v, bv0, offset_v) + call basis1(kv, kv+mv, x(2), knot_v, bv1, offset_v) + call basis2(kv, kv+mv, x(2), knot_v, bv2, offset_v) + + P(:) = 0.0 + Pu(:) = 0.0 + Pv(:) = 0.0 + Puu(:) = 0.0 + Puv(:) = 0.0 + Pvv(:) = 0.0 + + do lu = 1, ku + hu = lu + offset_u + do lv = 1, kv + hv = lv + offset_v + P = P + bu0(lu) * bv0(lv) * cp(hu, hv, :) + Pu = Pu + bu1(lu) * bv0(lv) * cp(hu, hv, :) + Pv = Pv + bu0(lu) * bv1(lv) * cp(hu, hv, :) + Puu = Puu + bu2(lu) * bv0(lv) * cp(hu, hv, :) + Puv = Puv + bu1(lu) * bv1(lv) * cp(hu, hv, :) + Pvv = Pvv + bu0(lu) * bv2(lv) * cp(hu, hv, :) + end do + end do + + f = P - pt0 + g(1) = 2 * dot_product(f, Pu) + g(2) = 2 * dot_product(f, Pv) + H(1, 1) = 2 * dot_product(Pu, Pu) + 2 * dot_product(f, Puu) + H(1, 2) = 2 * dot_product(Pu, Pv) + 2 * dot_product(f, Puv) + H(2, 2) = 2 * dot_product(Pv, Pv) + 2 * dot_product(f, Pvv) + H(2, 1) = H(1, 2) + do k = 1, 2 + if (((x(k) .eq. 0) .and. (g(k) .gt. 0)) .or. & + ((x(k) .eq. 1) .and. (g(k) .lt. 0))) then + g(k) = 0.0 + H(1, 2) = 0.0 + H(2, 1) = 0.0 + H(k, k) = 1.0 + end if + end do + det = H(1, 1) * H(2, 2) - H(1, 2) * H(2, 1) + W(1, 1) = H(2, 2) / det + W(1, 2) = -H(1, 2) / det + W(2, 2) = H(1, 1) / det + W(2, 1) = W(1, 2) + norm_g = sqrt(dot_product(g, g)) + dx(1) = -dot_product(W(1, :), g) + dx(2) = -dot_product(W(2, :), g) + do k = 1, 2 + if (x(k) + dx(k) .lt. 0) then + dx(k) = -x(k) + else if (x(k) + dx(k) .gt. 1) then + dx(k) = 1 - x(k) + end if + end do + norm_dx = sqrt(dot_product(dx, dx)) + norm_f = sqrt(dot_product(f, f)) + + ! print *, counter, norm + if (((norm_g .lt. 1e-13) .or. (norm_dx .lt. 1e-13)) .and. & + (norm_f .lt. 1e-6)) then + fail = .False. + exit + end if + x = x + dx + + end do + pt(:) = P(:) + +end subroutine newtonprojection + + + +subroutine cross_product(veca, vecb, vecc) + + implicit none + + double precision, intent(in) :: veca(3), vecb(3) + double precision, intent(out) :: vecc(3) + + vecc(3) = veca(1) * vecb(2) - veca(2) * vecb(1) + vecc(2) = veca(3) * vecb(1) - veca(1) * vecb(3) + vecc(1) = veca(2) * vecb(3) - veca(3) * vecb(2) + +end subroutine cross_product diff --git a/src/BSE/compute_pt_jacobian.f90 b/src/BSE/compute_pt_jacobian.f90 old mode 100644 new mode 100755 diff --git a/src/BSE/compute_sc_jacobian.f90 b/src/BSE/compute_sc_jacobian.f90 old mode 100644 new mode 100755 diff --git a/src/BSE/compute_topology.f90 b/src/BSE/compute_topology.f90 old mode 100644 new mode 100755 diff --git a/src/PGM/interpolant/computeCone.f90 b/src/PGM/interpolant/computeCone.f90 old mode 100644 new mode 100755 diff --git a/src/PGM/interpolant/computeJunction.f90 b/src/PGM/interpolant/computeJunction.f90 old mode 100644 new mode 100755 diff --git a/src/PGM/interpolant/computeTip.f90 b/src/PGM/interpolant/computeTip.f90 old mode 100644 new mode 100755 diff --git a/src/PGM/interpolant/interpolant.f90 b/src/PGM/interpolant/interpolant.f90 old mode 100644 new mode 100755 diff --git a/src/PGM/parameter/bspline_basis.f90 b/src/PGM/parameter/bspline_basis.f90 old mode 100644 new mode 100755 diff --git a/src/PGM/parameter/bspline_knot.f90 b/src/PGM/parameter/bspline_knot.f90 old mode 100644 new mode 100755 diff --git a/src/PGM/parameter/bspline_param.f90 b/src/PGM/parameter/bspline_param.f90 old mode 100644 new mode 100755 diff --git a/src/PGM/parameter/compute_bspline.f90 b/src/PGM/parameter/compute_bspline.f90 old mode 100644 new mode 100755 diff --git a/src/PGM/primitive/computeAngles.f90 b/src/PGM/primitive/computeAngles.f90 old mode 100644 new mode 100755 diff --git a/src/PGM/primitive/computeRotations.f90 b/src/PGM/primitive/computeRotations.f90 old mode 100644 new mode 100755 diff --git a/src/PGM/primitive/computeRtnMtx.f90 b/src/PGM/primitive/computeRtnMtx.f90 old mode 100644 new mode 100755 diff --git a/src/PGM/primitive/computeSections.f90 b/src/PGM/primitive/computeSections.f90 old mode 100644 new mode 100755 diff --git a/src/PGM/primitive/computeShape.f90 b/src/PGM/primitive/computeShape.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/BLS/assembleMtx.f90 b/src/PSM/BLS/assembleMtx.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/CDT/addNode.f90 b/src/PSM/CDT/addNode.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/CDT/computeCDT.f90 b/src/PSM/CDT/computeCDT.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/CDT/constraints.f90 b/src/PSM/CDT/constraints.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/CDT/delaunay.f90 b/src/PSM/CDT/delaunay.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/CDT/delete.f90 b/src/PSM/CDT/delete.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/CDT/misc.f90 b/src/PSM/CDT/misc.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/CDT/nearest.f90 b/src/PSM/CDT/nearest.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/CDT/output.f90 b/src/PSM/CDT/output.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/CDT/postProcess.f90 b/src/PSM/CDT/postProcess.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/computeAdjoiningEdges.f90 b/src/PSM/GFEM/computeAdjoiningEdges.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/computeEdgeLengths.f90 b/src/PSM/GFEM/computeEdgeLengths.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/computeFaceDimensions.f90 b/src/PSM/GFEM/computeFaceDimensions.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/computeFaceEdges.f90 b/src/PSM/GFEM/computeFaceEdges.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/computeFaceIntersections.f90 b/src/PSM/GFEM/computeFaceIntersections.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/computeGroupIntersections.f90 b/src/PSM/GFEM/computeGroupIntersections.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/computeGroupSplits.f90 b/src/PSM/GFEM/computeGroupSplits.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/computeIntersectionVerts.f90 b/src/PSM/GFEM/computeIntersectionVerts.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/computeMemberEdges.f90 b/src/PSM/GFEM/computeMemberEdges.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/computeMemberNodes.f90 b/src/PSM/GFEM/computeMemberNodes.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/computeMemberTopology.f90 b/src/PSM/GFEM/computeMemberTopology.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/computeMembers.f90 b/src/PSM/GFEM/computeMembers.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/computePreviewMembers.f90 b/src/PSM/GFEM/computePreviewMembers.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/computePreviewSurfaces.f90 b/src/PSM/GFEM/computePreviewSurfaces.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/computeProjtnInputs.f90 b/src/PSM/GFEM/computeProjtnInputs.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/computeSurfaceProjections.f90 b/src/PSM/GFEM/computeSurfaceProjections.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/computeSurfaces.f90 b/src/PSM/GFEM/computeSurfaces.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/identifySymmNodes.f90 b/src/PSM/GFEM/identifySymmNodes.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/importMembers.f90 b/src/PSM/GFEM/importMembers.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/misc.f90 b/src/PSM/GFEM/misc.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/removeDuplicateNodes.f90 b/src/PSM/GFEM/removeDuplicateNodes.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/GFEM/removeRightQuads.f90 b/src/PSM/GFEM/removeRightQuads.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/QUAD/addEdgePts.f90 b/src/PSM/QUAD/addEdgePts.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/QUAD/addInteriorPts.f90 b/src/PSM/QUAD/addInteriorPts.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/QUAD/addIntersectionPts.f90 b/src/PSM/QUAD/addIntersectionPts.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/QUAD/computeAdjMap.f90 b/src/PSM/QUAD/computeAdjMap.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/QUAD/computeConstraints.f90 b/src/PSM/QUAD/computeConstraints.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/QUAD/computeDivisions.f90 b/src/PSM/QUAD/computeDivisions.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/QUAD/computeGrid.f90 b/src/PSM/QUAD/computeGrid.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/QUAD/computeQuad2Edge.f90 b/src/PSM/QUAD/computeQuad2Edge.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/QUAD/computeQuadDominant.f90 b/src/PSM/QUAD/computeQuadDominant.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/QUAD/computeQuads.f90 b/src/PSM/QUAD/computeQuads.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/QUAD/computeTriangles.f90 b/src/PSM/QUAD/computeTriangles.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/QUAD/importEdges.f90 b/src/PSM/QUAD/importEdges.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/QUAD/removeDegenerateEdges.f90 b/src/PSM/QUAD/removeDegenerateEdges.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/QUAD/removeDuplicateEdges.f90 b/src/PSM/QUAD/removeDuplicateEdges.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/QUAD/removeDuplicateQuads.f90 b/src/PSM/QUAD/removeDuplicateQuads.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/QUAD/removeDuplicateTriangles.f90 b/src/PSM/QUAD/removeDuplicateTriangles.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/QUAD/removeDuplicateVerts.f90 b/src/PSM/QUAD/removeDuplicateVerts.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/QUAD/removeInvalidQuads.f90 b/src/PSM/QUAD/removeInvalidQuads.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/QUAD/reorderCollinear.f90 b/src/PSM/QUAD/reorderCollinear.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/QUAD/splitEdges.f90 b/src/PSM/QUAD/splitEdges.f90 old mode 100644 new mode 100755 diff --git a/src/PSM/QUAD/splitTrisNQuads.f90 b/src/PSM/QUAD/splitTrisNQuads.f90 old mode 100644 new mode 100755