From 515b7128237ccaf11abbb90544c1b85dc4beb10f Mon Sep 17 00:00:00 2001 From: Mike Sprague Date: Tue, 6 Oct 2020 15:40:21 -0600 Subject: [PATCH 1/8] Added test for BD_diffmtc (5-node problem) Fixed call to diffmtc in original 2-node problem --- modules/beamdyn/tests/test_BD_diffmtc.F90 | 213 +++++++++++++++------- 1 file changed, 149 insertions(+), 64 deletions(-) diff --git a/modules/beamdyn/tests/test_BD_diffmtc.F90 b/modules/beamdyn/tests/test_BD_diffmtc.F90 index 7b4bc91889..ca5810d26a 100644 --- a/modules/beamdyn/tests/test_BD_diffmtc.F90 +++ b/modules/beamdyn/tests/test_BD_diffmtc.F90 @@ -1,15 +1,12 @@ -@test -subroutine test_BD_diffmtc() - ! branches to test - ! - 2 node, 1 element - +module test_BD_diffmtc + use pFUnit_mod use BeamDyn use NWTC_Num use test_tools - + implicit none - + integer :: n, i type(BD_ParameterType) :: parametertype real(BDKi), allocatable :: test_shape(:,:), test_shapederivative(:,:) @@ -19,60 +16,148 @@ subroutine test_BD_diffmtc() character :: ErrMsg character(1024) :: testname real(BDKi) :: tolerance - - ! initialize NWTC_Num constants - call SetConstants() - - tolerance = 1e-14 - - - ! -------------------------------------------------------------------------- - testname = "2 node, 1 element; all quadrature points are at GLL nodes:" - - ! the shape functions should be: - ! h1(-1) = 1, h1(+1) = 0 - ! h2(-1) = 0, h2(+1) = 1 - ! - ! this is satisfied by these linear equations - ! h1(s) = 0.5*(1-s) - ! h2(s) = 0.5*(1+s) - ! therefore, - ! h1` = -0.5 - ! h2` = 0.5 - ! - ! the expected result of BD_diffmtc is - ! Shp - the shape function evaluated at the GLL nodes - ! ShpDer - the shape function derivative evaluated at the GLL nodes - ! - ! shp(1,:) = 1.0, 0.0 - ! shp(2,:) = 0.0, 1.0 - ! shpder(1,:) = -0.5, -0.5 - ! shpder(2,:) = 0.5, 0.5 - - parametertype = simpleParameterType() - parametertype%nodes_per_elem = 2 - parametertype%nqp = 2 - n = parametertype%nodes_per_elem - - call AllocAry(test_shape, parametertype%nodes_per_elem, parametertype%nqp, "test_shape", ErrStat, ErrMsg) - call AllocAry(test_shapederivative, parametertype%nodes_per_elem, parametertype%nqp, "test_shapederivative", ErrStat, ErrMsg) - - call AllocAry(parametertype%QPtN, parametertype%nodes_per_elem, 'QPtN', ErrStat, ErrMsg) - parametertype%QPtN = (/ -1.0, 1.0 /) - - call AllocAry(gll_nodes, n, "GLL points array", ErrStat, ErrMsg) - gll_nodes = (/ -1.0, 1.0 /) - - call BD_diffmtc(parametertype%nodes_per_elem, gll_nodes, parametertype%QPtN, parametertype%nqp, test_shape, test_shapederivative) - - call AllocAry(baseline_shape, parametertype%nqp, parametertype%nodes_per_elem, "baseline_shape", ErrStat, ErrMsg) - call AllocAry(baseline_shapederivative, parametertype%nqp, parametertype%nodes_per_elem, "baseline_shapederivative", ErrStat, ErrMsg) - baseline_shape(1,:) = (/ 1.0, 0.0 /) - baseline_shape(2,:) = (/ 0.0, 1.0 /) - baseline_shapederivative(1,:) = (/ -0.5, -0.5 /) - baseline_shapederivative(2,:) = (/ 0.5, 0.5 /) - - @assertEqual(baseline_shape, test_shape, tolerance, testname) - @assertEqual(baseline_shapederivative, test_shapederivative, tolerance, testname) - -end subroutine + + ! mathematica code for calculating shape function values and derivative values + !ClearAll[p, \[Xi]j, xj, h, j, \[Xi], x, sol]; + !(* find p+1 GLL points *) + !p = 4; + !test[x_] = (1 - x^2) D[LegendreP[p, x], x]; + !\[Xi]j = Array[f, p + 1]; (* initialize list *) + !For[j = 0, j <= p, j++, + !sol = x /. FindRoot[test[x], {x, -Cos[(Pi (j))/p]}]; + ! \[Xi]j[[j + 1]] = sol + !] + !h[\[Xi]_, xj_] = -(((1. - \[Xi]^2) D[ LegendreP[p, \[Xi]], \[Xi]])/( + !p (p + 1) LegendreP[p, xj] (\[Xi] - xj))); + !hd[\[Xi]_, xj_] = D[h[\[Xi], xj], \[Xi]] + !Plot[h[x, \[Xi]j[[1]]], {x, -1, 1}, PlotRange -> All] + !h[0.5, \[Xi]j[[2]]] + +contains + + @test + subroutine test_BD_diffmtc_2node() + ! branches to test + ! - 2 node, 1 element + + ! initialize NWTC_Num constants + call SetConstants() + + tolerance = 1e-14 + + ! -------------------------------------------------------------------------- + testname = "2-node element: evaluate shape/shapederivative at nodes" + + ! the shape functions should be: + ! h1(-1) = 1, h1(+1) = 0 + ! h2(-1) = 0, h2(+1) = 1 + ! + ! this is satisfied by these linear equations + ! h1(s) = 0.5*(1-s) + ! h2(s) = 0.5*(1+s) + ! therefore, + ! h1` = -0.5 + ! h2` = 0.5 + ! + ! the expected result of BD_diffmtc is + ! Shp - the shape function evaluated at the GLL nodes + ! ShpDer - the shape function derivative evaluated at the GLL nodes + ! + ! shp(1,:) = 1.0, 0.0 + ! shp(2,:) = 0.0, 1.0 + ! shpder(1,:) = -0.5, -0.5 + ! shpder(2,:) = 0.5, 0.5 + + parametertype = simpleParameterType() + parametertype%nodes_per_elem = 2 + parametertype%nqp = 2 + n = parametertype%nodes_per_elem + + call AllocAry(test_shape, parametertype%nodes_per_elem, parametertype%nqp, "test_shape", ErrStat, ErrMsg) + call AllocAry(test_shapederivative, parametertype%nodes_per_elem, parametertype%nqp, "test_shapederivative", ErrStat, ErrMsg) + + call AllocAry(parametertype%QPtN, parametertype%nodes_per_elem, 'QPtN', ErrStat, ErrMsg) + parametertype%QPtN = (/ -1.0, 1.0 /) + + call AllocAry(gll_nodes, n, "GLL points array", ErrStat, ErrMsg) + gll_nodes = (/ -1.0, 1.0 /) + + call BD_diffmtc(parametertype%nodes_per_elem, gll_nodes, parametertype%QPtN, parametertype%nqp, test_shape, test_shapederivative) + + call AllocAry(baseline_shape, parametertype%nodes_per_elem,parametertype%nqp, "baseline_shape", ErrStat, ErrMsg) + call AllocAry(baseline_shapederivative, parametertype%nodes_per_elem,parametertype%nqp, "baseline_shapederivative", ErrStat, ErrMsg) + baseline_shape(1,:) = (/ 1.0, 0.0 /) + baseline_shape(2,:) = (/ 0.0, 1.0 /) + baseline_shapederivative(1,:) = (/ -0.5, -0.5 /) + baseline_shapederivative(2,:) = (/ 0.5, 0.5 /) + + @assertEqual(baseline_shape, test_shape, tolerance, testname) + @assertEqual(baseline_shapederivative, test_shapederivative, tolerance, testname) + + deallocate(test_shape) + deallocate(test_shapederivative) + deallocate(parametertype%QPtN) + deallocate(gll_nodes) + deallocate(baseline_shape) + deallocate(baseline_shapederivative) + + end subroutine + + @test + subroutine test_BD_diffmtc_5node() + ! branches to test + ! - 5 node, 1 element + + ! initialize NWTC_Num constants + call SetConstants() + + tolerance = 1e-14 + + ! -------------------------------------------------------------------------- + testname = "5-node element: evaluate shape/shapederivative at nodes and at three non-node locations" + + parametertype = simpleParameterType() + parametertype%nodes_per_elem = 5 + parametertype%nqp = 8 ! testing the GLL nodes and three non-nodal locations (-0.8, 0.1, 0.4) + n = parametertype%nodes_per_elem + + call AllocAry(test_shape, parametertype%nodes_per_elem, parametertype%nqp, "test_shape", ErrStat, ErrMsg) + call AllocAry(test_shapederivative, parametertype%nodes_per_elem, parametertype%nqp, "test_shapederivative", ErrStat, ErrMsg) + + call AllocAry(parametertype%QPtN, parametertype%nodes_per_elem, 'QPtN', ErrStat, ErrMsg) + ! in following, first five points are the 4th-order GLL locations; last three points are randomly chosen off-node points + parametertype%QPtN = (/ -1.0, -0.6546536707079771, 0.0, 0.6546536707079771, 1.0, -0.8, 0.1, 0.4 /) + + call AllocAry(gll_nodes, n, "GLL points array", ErrStat, ErrMsg) + gll_nodes = (/ -1.0, -0.6546536707079771, 0.0, 0.6546536707079771, 1.0 /) + + call BD_diffmtc(parametertype%nodes_per_elem, gll_nodes, parametertype%QPtN, parametertype%nqp, test_shape, test_shapederivative) + + call AllocAry(baseline_shape, parametertype%nodes_per_elem, parametertype%nqp, "baseline_shape", ErrStat, ErrMsg) + call AllocAry(baseline_shapederivative, parametertype%nodes_per_elem, parametertype%nqp, "baseline_shapederivative", ErrStat, ErrMsg) + + ! baseline_shape and baseline_shapederivative were calculated in Mathematica for fourth order Legendre spectral FE + + baseline_shape(1,:) = (/ 1., 0., 0., 0., 0., 0.266400000000000, 0.0329625, 0.0564/) + baseline_shape(2,:) = (/ 0., 1., 0., 0., 0., 0.855336358376291, -0.1121093731918499, -0.1746924181056723/) + baseline_shape(3,:) = (/ 0., 0., 1., 0., 0., -0.1776000000000001, 0.9669, 0.5263999999999998/) + baseline_shape(4,:) = (/ 0., 0., 0., 1., 0., 0.0854636416237095, 0.1525343731918499, 0.7234924181056726/) + baseline_shape(5,:) = (/ 0., 0., 0., 0., 1., -0.0296000000000000, -0.0402875, -0.1316/) + baseline_shapederivative(1,:) = (/-5., -1.240990253030983,0.375,-0.2590097469690172,0.5,-2.497, 0.27725, -0.1210000000000001/) + baseline_shapederivative(2,:) = (/ 6.756502488724241,0.,-1.336584577695454,0.7637626158259736,-1.410164177942427,2.144324478146486,-0.896320373697923, 0.4156426862650312/) + baseline_shapederivative(3,:) = (/-2.666666666666667, 1.74574312188794,0.,-1.74574312188794,2.666666666666667,0.5546666666666656,-0.6573333333333338, -2.069333333333333/) + baseline_shapederivative(4,:) = (/ 1.410164177942427,-0.7637626158259736, 1.336584577695454,0.,-6.756502488724241,-0.31499114481315,1.696653707031257, 1.805690647068303/) + baseline_shapederivative(5,:) = (/-0.5, 0.2590097469690172,-0.375,1.240990253030983,5.,0.1129999999999999,-0.42025, -0.03099999999999978/) + + @assertEqual(baseline_shape, test_shape, tolerance, testname) + @assertEqual(baseline_shapederivative, test_shapederivative, tolerance, testname) + + deallocate(test_shape) + deallocate(test_shapederivative) + deallocate(parametertype%QPtN) + deallocate(gll_nodes) + deallocate(baseline_shape) + deallocate(baseline_shapederivative) + + end subroutine +end module From a240c01bf7ee933bef73634cdbe91fa93d9201a9 Mon Sep 17 00:00:00 2001 From: Mike Sprague Date: Thu, 8 Oct 2020 13:07:57 -0600 Subject: [PATCH 2/8] pushing test_BD_InitShpDerJaco.F90 for help from raf --- .../beamdyn/tests/test_BD_InitShpDerJaco.F90 | 478 ++++++++++++------ 1 file changed, 326 insertions(+), 152 deletions(-) diff --git a/modules/beamdyn/tests/test_BD_InitShpDerJaco.F90 b/modules/beamdyn/tests/test_BD_InitShpDerJaco.F90 index 184ca3ec1b..cfba996cbf 100644 --- a/modules/beamdyn/tests/test_BD_InitShpDerJaco.F90 +++ b/modules/beamdyn/tests/test_BD_InitShpDerJaco.F90 @@ -1,8 +1,5 @@ -@test -subroutine test_BD_InitShpDerJaco() - ! branches to test - ! - 3 node, 1 element; undeformed - +module test_BD_InitShpDerJaco + use pFUnit_mod use BeamDyn use NWTC_Num @@ -13,6 +10,7 @@ subroutine test_BD_InitShpDerJaco() integer(IntKi) :: i, j, idx_qp, nelem type(BD_ParameterType) :: p real(BDKi), allocatable :: gll_nodes(:), inp_QPtWeight(:) + real(BDKi), allocatable :: baseline_QPtWeight(:), baseline_QPtN(:) real(BDKi), allocatable :: baseline_Shp(:,:), baseline_ShpDer(:,:), baseline_jacobian(:,:), baseline_QPtw_ShpDer(:,:) real(BDKi), allocatable :: baseline_QPtw_Shp_ShpDer(:,:,:), baseline_QPtw_Shp_Jac(:,:,:) real(BDKi), allocatable :: baseline_QPtw_Shp_Shp_Jac(:,:,:,:), baseline_QPtw_ShpDer_ShpDer_Jac(:,:,:,:) @@ -21,171 +19,347 @@ subroutine test_BD_InitShpDerJaco() character(1024) :: testname real(BDKi) :: tolerance - ! initialize NWTC_Num constants - call SetConstants() - - tolerance = 1e-14 - - - ! -------------------------------------------------------------------------- - testname = "3 node, 1 element, undeformed:" - - ! Lets assume a 3 node element with parametric coordinate s where s lies in [-1 1] - ! The three nodes lie at s1 = -1, s2 = 0, and s3 = 1 - ! The corresponding positions of the nodes in physical space are assumed to be - ! x1 = 0, x2 = 10, and x2 = 20 - ! - ! The quadratic shape functions for such a 3 node element at any given 's' are given by - ! N1(s) = (s2-s)*(s3-s)/(s2-s1)/(s3-s1) - ! N2(s) = (s1-s)*(s3-s)/(s1-s2)/(s3-s2) - ! N3(s) = (s1-s)*(s2-s)/(s1-s3)/(s2-s3) - ! - ! We are interested in testing the function test_BD_InitShpDerJaco for s = 0.5 - ! At s = 0.5, the above values are computed to be N1 = -0.125, N2 = 0.75, and N3 = -0.375 - ! - ! The derivatives dN/ds are computed to be dN1 = -, dN2 = -1, and dN3 = 1 - ! - ! the Jacobian is invariant of the quadrature point and account for change in volume - ! between physical and parametric space. Jacobian = meas(x)/meas(s), where meas is the - ! measured quantity (volume in 3D, area in 2D, and lenght in 1D) - ! The Jacobian for this element is given by 20/2 = 10 - - ! build the p object based on the above mentioned test model - p = simpleparametertype() - p%elem_total = 1 - p%nodes_per_elem = 3 - p%nqp = 1 - - ! Allocate memory for baseline results - call AllocAry(baseline_Shp , p%nodes_per_elem, p%nqp, 'Reference Shp' , ErrStat, ErrMsg) - call AllocAry(baseline_ShpDer , p%nodes_per_elem, p%nqp, 'Reference ShpDer' , ErrStat, ErrMsg) - call AllocAry(baseline_Jacobian, p%elem_total , p%nqp, 'Reference Jacobian', ErrStat, ErrMsg) - call AllocAry(inp_QPtWeight , p%nqp , 'QPtWeight' , ErrStat, ErrMsg) - - ! Allocate memory for other relevant variables belonging to module p - call AllocAry(baseline_QPtw_Shp_Shp_Jac , p%nqp, p%nodes_per_elem, p%nodes_per_elem, p%elem_total, 'reference QPtw_Shp_Shp_Jac' , ErrStat, ErrMsg) - call AllocAry(baseline_QPtw_ShpDer_ShpDer_Jac, p%nqp, p%nodes_per_elem, p%nodes_per_elem, p%elem_total, 'reference baseline_QPtw_ShpDer_ShpDer_Jac', ErrStat, ErrMsg) - call AllocAry(baseline_QPtw_Shp_ShpDer , p%nqp, p%nodes_per_elem, p%nodes_per_elem , 'reference QPtw_Shp_ShpDer' , ErrStat, ErrMsg) - call AllocAry(baseline_QPtw_Shp_Jac , p%nqp, p%nodes_per_elem, p%elem_total , 'reference QPtw_Shp_Jac' , ErrStat, ErrMsg) - call AllocAry(baseline_QPtw_ShpDer , p%nqp, p%nodes_per_elem , 'reference QPtw_ShpDer' , ErrStat, ErrMsg) - - ! assign baseling results - ! assign baseline jacobian based on example as described above - baseline_jacobian(1,1) = 10.0 ! we assume 1 element quadrature point. Hence we set only the index (1,1) - - ! assign baseline shape functions based on example as described above - baseline_Shp(1,1) = -0.125 - baseline_Shp(2,1) = 0.750 - baseline_Shp(3,1) = 0.375 - - ! assign baseline shape function derivatives based on example as described above - baseline_ShpDer(1,1) = 0.0 - baseline_ShpDer(2,1) = -1.0 - baseline_ShpDer(3,1) = 1.0 +contains - ! assign the weight to the quadrature point which is an input parameter - inp_QPtWeight(1) = 1.0; - - ! Allocate memory for relevant variables belonging to module p - call AllocAry(p%Shp , p%nodes_per_elem, p%nqp, 'Shp' , ErrStat, ErrMsg) - call AllocAry(p%ShpDer , p%nodes_per_elem, p%nqp, 'ShpDer' , ErrStat, ErrMsg) - call AllocAry(p%uuN0, 3, p%nodes_per_elem, p%nqp, 'uuN0' , ErrStat, ErrMsg) - call AllocAry(p%Jacobian , p%elem_total , p%nqp, 'Jacobian' , ErrStat, ErrMsg) - call AllocAry(p%QPtN , p%nodes_per_elem , 'QPtN' , ErrStat, ErrMsg) - call AllocAry(p%QPtWeight, p%nqp , 'QPtWeight', ErrStat, ErrMsg) - - ! Allocate memory for other relevant variables belonging to module p - call AllocAry(p%QPtw_Shp_Shp_Jac , p%nqp, p%nodes_per_elem, p%nodes_per_elem,p%elem_total, 'QPtw_Shp_Shp_Jac', ErrStat, ErrMsg) - call AllocAry(p%QPtw_Shp_ShpDer , p%nqp, p%nodes_per_elem, p%nodes_per_elem , 'QPtw_Shp_ShpDer', ErrStat, ErrMsg) - call AllocAry(p%QPtw_ShpDer_ShpDer_Jac, p%nqp, p%nodes_per_elem, p%nodes_per_elem,p%elem_total, 'QPtw_ShpDer_ShpDer_Jac', ErrStat, ErrMsg) - call AllocAry(p%QPtw_Shp_Jac , p%nqp, p%nodes_per_elem, p%elem_total , 'QPtw_Shp_Jac', ErrStat, ErrMsg) - call AllocAry(p%QPtw_ShpDer , p%nqp, p%nodes_per_elem , 'QPtw_ShpDer', ErrStat, ErrMsg) - - ! uuN0 is of dimension (3 dof, nodes_per_elem, elem_total) - p%uuN0(1:3,1,1) = (/ 0.0, 0.0, 0.0 /) - p%uuN0(1:3,2,1) = (/ 0.0, 0.0, 10.0 /) - p%uuN0(1:3,3,1) = (/ 0.0, 0.0, 20.0 /) + @test + subroutine test_BD_InitShpDerJaco_3node() + ! branches to test + ! - 3 node, 1 element; undeformed + - p%QPtN = (/ 0.5 /) ! Note, we assume 1 quadrature point + ! initialize NWTC_Num constants + call SetConstants() + + tolerance = 1e-14 + + + ! -------------------------------------------------------------------------- + testname = "3 node, 1 element, undeformed:" + + ! Lets assume a 3 node element with parametric coordinate s where s lies in [-1 1] + ! The three nodes lie at s1 = -1, s2 = 0, and s3 = 1 + ! The corresponding positions of the nodes in physical space are assumed to be + ! x1 = 0, x2 = 10, and x2 = 20 + ! + ! The quadratic shape functions for such a 3 node element at any given 's' are given by + ! N1(s) = (s2-s)*(s3-s)/(s2-s1)/(s3-s1) + ! N2(s) = (s1-s)*(s3-s)/(s1-s2)/(s3-s2) + ! N3(s) = (s1-s)*(s2-s)/(s1-s3)/(s2-s3) + ! + ! We are interested in testing the function test_BD_InitShpDerJaco for s = 0.5 + ! At s = 0.5, the above values are computed to be N1 = -0.125, N2 = 0.75, and N3 = -0.375 + ! + ! The derivatives dN/ds are computed to be dN1 = -, dN2 = -1, and dN3 = 1 + ! + ! the Jacobian is invariant of the quadrature point and account for change in volume + ! between physical and parametric space. Jacobian = meas(x)/meas(s), where meas is the + ! measured quantity (volume in 3D, area in 2D, and lenght in 1D) + ! The Jacobian for this element is given by 20/2 = 10 - p%QPtWeight = inp_QPtWeight ! Since we assume 1 quadrature point, the weight by defualt = 1 - - ! Allocate memory for GLL node positions in 1D parametric space - call AllocAry(gll_nodes, p%nodes_per_elem, "GLL points array", ErrStat, ErrMsg) - gll_nodes = (/ -1.0, 0.0, 1.0 /) - - ! call the test subroutine - call BD_InitShpDerJaco(gll_nodes, p) - - ! check the baseline shape functions and their derivatives - do idx_qp = 1, p%nqp - do j = 1, p%nodes_per_elem - @assertEqual(baseline_Shp(j,idx_qp) , p%Shp(j,idx_qp) , tolerance, testname) - @assertEqual(baseline_ShpDer(j,idx_qp), p%ShpDer(j,idx_qp), tolerance, testname) - end do - end do - - ! check the baseline jacobian - do nelem = 1, p%elem_total + ! build the p object based on the above mentioned test model + p = simpleparametertype() + p%elem_total = 1 + p%nodes_per_elem = 3 + p%nqp = 1 + + ! Allocate memory for baseline results FIXME -- seems like arguments reversed + call AllocAry(baseline_Shp , p%nodes_per_elem, p%nqp, 'Reference Shp' , ErrStat, ErrMsg) + call AllocAry(baseline_ShpDer , p%nodes_per_elem, p%nqp, 'Reference ShpDer' , ErrStat, ErrMsg) + call AllocAry(baseline_Jacobian, p%elem_total , p%nqp, 'Reference Jacobian', ErrStat, ErrMsg) + call AllocAry(inp_QPtWeight , p%nqp , 'QPtWeight' , ErrStat, ErrMsg) + + ! Allocate memory for other relevant variables belonging to module p + call AllocAry(baseline_QPtw_Shp_Shp_Jac , p%nqp, p%nodes_per_elem, p%nodes_per_elem, p%elem_total, 'reference QPtw_Shp_Shp_Jac' , ErrStat, ErrMsg) + call AllocAry(baseline_QPtw_ShpDer_ShpDer_Jac, p%nqp, p%nodes_per_elem, p%nodes_per_elem, p%elem_total, 'reference baseline_QPtw_ShpDer_ShpDer_Jac', ErrStat, ErrMsg) + call AllocAry(baseline_QPtw_Shp_ShpDer , p%nqp, p%nodes_per_elem, p%nodes_per_elem , 'reference QPtw_Shp_ShpDer' , ErrStat, ErrMsg) + call AllocAry(baseline_QPtw_Shp_Jac , p%nqp, p%nodes_per_elem, p%elem_total , 'reference QPtw_Shp_Jac' , ErrStat, ErrMsg) + call AllocAry(baseline_QPtw_ShpDer , p%nqp, p%nodes_per_elem , 'reference QPtw_ShpDer' , ErrStat, ErrMsg) + + ! assign baseling results + ! assign baseline jacobian based on example as described above + baseline_jacobian(1,1) = 10.0 ! we assume 1 element quadrature point. Hence we set only the index (1,1) + + ! assign baseline shape functions based on example as described above + baseline_Shp(1,1) = -0.125 + baseline_Shp(2,1) = 0.750 + baseline_Shp(3,1) = 0.375 + + ! assign baseline shape function derivatives based on example as described above + baseline_ShpDer(1,1) = 0.0 + baseline_ShpDer(2,1) = -1.0 + baseline_ShpDer(3,1) = 1.0 + + ! assign the weight to the quadrature point which is an input parameter + inp_QPtWeight(1) = 1.0; + + ! Allocate memory for relevant variables belonging to module p FIXME -- seems like arguments reversed; work becuase nodes_per_elem = nqp = 1 + call AllocAry(p%Shp , p%nodes_per_elem, p%nqp, 'Shp' , ErrStat, ErrMsg) + call AllocAry(p%ShpDer , p%nodes_per_elem, p%nqp, 'ShpDer' , ErrStat, ErrMsg) + call AllocAry(p%uuN0, 3, p%nodes_per_elem, p%elem_total, 'uuN0' , ErrStat, ErrMsg) + call AllocAry(p%Jacobian , p%elem_total , p%nqp, 'Jacobian' , ErrStat, ErrMsg) ! reversed arguments? + call AllocAry(p%QPtN , p%nodes_per_elem , 'QPtN' , ErrStat, ErrMsg) + call AllocAry(p%QPtWeight, p%nqp , 'QPtWeight', ErrStat, ErrMsg) + + ! Allocate memory for other relevant variables belonging to module p + call AllocAry(p%QPtw_Shp_Shp_Jac , p%nqp, p%nodes_per_elem, p%nodes_per_elem,p%elem_total, 'QPtw_Shp_Shp_Jac', ErrStat, ErrMsg) + call AllocAry(p%QPtw_Shp_ShpDer , p%nqp, p%nodes_per_elem, p%nodes_per_elem , 'QPtw_Shp_ShpDer', ErrStat, ErrMsg) + call AllocAry(p%QPtw_ShpDer_ShpDer_Jac, p%nqp, p%nodes_per_elem, p%nodes_per_elem,p%elem_total, 'QPtw_ShpDer_ShpDer_Jac', ErrStat, ErrMsg) + call AllocAry(p%QPtw_Shp_Jac , p%nqp, p%nodes_per_elem, p%elem_total , 'QPtw_Shp_Jac', ErrStat, ErrMsg) + call AllocAry(p%QPtw_ShpDer , p%nqp, p%nodes_per_elem , 'QPtw_ShpDer', ErrStat, ErrMsg) + + ! uuN0 is of dimension (3 dof, nodes_per_elem, elem_total) + p%uuN0(1:3,1,1) = (/ 0.0, 0.0, 0.0 /) + p%uuN0(1:3,2,1) = (/ 0.0, 0.0, 10.0 /) + p%uuN0(1:3,3,1) = (/ 0.0, 0.0, 20.0 /) + + p%QPtN = (/ 0.5 /) ! Note, we assume 1 quadrature point + + p%QPtWeight = inp_QPtWeight ! Since we assume 1 quadrature point, the weight by defualt = 1 + + ! Allocate memory for GLL node positions in 1D parametric space + call AllocAry(gll_nodes, p%nodes_per_elem, "GLL points array", ErrStat, ErrMsg) + gll_nodes = (/ -1.0, 0.0, 1.0 /) + + ! call the test subroutine + call BD_InitShpDerJaco(gll_nodes, p) + + ! check the baseline shape functions and their derivatives do idx_qp = 1, p%nqp - @assertEqual(baseline_jacobian(nelem,idx_qp), p%jacobian(nelem,idx_qp), tolerance, testname) + do j = 1, p%nodes_per_elem + @assertEqual(baseline_Shp(j,idx_qp) , p%Shp(j,idx_qp) , tolerance, testname) + @assertEqual(baseline_ShpDer(j,idx_qp), p%ShpDer(j,idx_qp), tolerance, testname) + end do end do - end do - - ! Test and assemble variables N*N^T*wt*Jacobian and dN*dN^T*wt/Jacobian - do nelem = 1, p%elem_total + + ! check the baseline jacobian + do nelem = 1, p%elem_total + do idx_qp = 1, p%nqp + @assertEqual(baseline_jacobian(nelem,idx_qp), p%jacobian(nelem,idx_qp), tolerance, testname) + end do + end do + + ! Test and assemble variables N*N^T*wt*Jacobian and dN*dN^T*wt/Jacobian + do nelem = 1, p%elem_total + do idx_qp = 1, p%nqp + do j = 1, p%nodes_per_elem + do i = 1, p%nodes_per_elem + ! Check the variable N*N^T*Jacobian + baseline_QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem) = baseline_Shp(i,idx_qp)*baseline_Shp(j,idx_qp)*inp_QPtWeight(idx_qp)*baseline_jacobian(idx_qp,nelem) + @assertEqual(baseline_QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem), p%QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem), tolerance, testname) + + ! Check the variable dN*dN^T*Jacobian + baseline_QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem) = baseline_ShpDer(i,idx_qp)*baseline_ShpDer(j,idx_qp)*inp_QPtWeight(idx_qp)/baseline_jacobian(idx_qp,nelem) + @assertEqual(baseline_QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem), p%QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem), tolerance, testname) + end do + end do + end do + end do + + ! Test and assemble variable N*dN^T*wt*Jacobian do idx_qp = 1, p%nqp do j = 1, p%nodes_per_elem do i = 1, p%nodes_per_elem - ! Check the variable N*N^T*Jacobian - baseline_QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem) = baseline_Shp(i,idx_qp)*baseline_Shp(j,idx_qp)*inp_QPtWeight(idx_qp)*baseline_jacobian(idx_qp,nelem) - @assertEqual(baseline_QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem), p%QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem), tolerance, testname) - - ! Check the variable dN*dN^T*Jacobian - baseline_QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem) = baseline_ShpDer(i,idx_qp)*baseline_ShpDer(j,idx_qp)*inp_QPtWeight(idx_qp)/baseline_jacobian(idx_qp,nelem) - @assertEqual(baseline_QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem), p%QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem), tolerance, testname) + baseline_QPtw_Shp_ShpDer(idx_qp,i,j) = baseline_Shp(i,idx_qp)*baseline_ShpDer(j,idx_qp)*inp_QPtWeight(idx_qp) + @assertEqual(baseline_QPtw_Shp_ShpDer(idx_qp,i,j), p%QPtw_Shp_ShpDer(idx_qp,i,j), tolerance, testname) end do end do end do - end do - - ! Test and assemble variable N*dN^T*wt*Jacobian - do idx_qp = 1, p%nqp - do j = 1, p%nodes_per_elem + + ! Test and assemble variable N*wt*Jacobian + do nelem = 1, p%elem_total do i = 1, p%nodes_per_elem - baseline_QPtw_Shp_ShpDer(idx_qp,i,j) = baseline_Shp(i,idx_qp)*baseline_ShpDer(j,idx_qp)*inp_QPtWeight(idx_qp) - @assertEqual(baseline_QPtw_Shp_ShpDer(idx_qp,i,j), p%QPtw_Shp_ShpDer(idx_qp,i,j), tolerance, testname) + do idx_qp = 1, p%nqp + baseline_QPtw_Shp_Jac(idx_qp,i,nelem) = baseline_Shp(i,idx_qp)*inp_QPtWeight(idx_qp)*baseline_Jacobian(idx_qp,nelem) + @assertEqual(baseline_QPtw_Shp_Jac(idx_qp,i,nelem), p%QPtw_Shp_Jac(idx_qp,i,nelem), tolerance, testname) + end do end do end do - end do - - ! Test and assemble variable N*wt*Jacobian - do nelem = 1, p%elem_total + + ! Test and assemble variable dN*wt. do i = 1, p%nodes_per_elem do idx_qp = 1, p%nqp - baseline_QPtw_Shp_Jac(idx_qp,i,nelem) = baseline_Shp(i,idx_qp)*inp_QPtWeight(idx_qp)*baseline_Jacobian(idx_qp,nelem) - @assertEqual(baseline_QPtw_Shp_Jac(idx_qp,i,nelem), p%QPtw_Shp_Jac(idx_qp,i,nelem), tolerance, testname) + baseline_QPtw_ShpDer(idx_qp,i) = baseline_ShpDer(i,idx_qp)*inp_QPtWeight(idx_qp) + @assertEqual(baseline_QPtw_ShpDer(idx_qp,i), p%QPtw_ShpDer(idx_qp,i), tolerance, testname) end do end do - end do + + ! dealocate baseline variables + deallocate(baseline_Shp) + deallocate(baseline_ShpDer) + deallocate(baseline_Jacobian) + deallocate(inp_QPtWeight) + deallocate(baseline_QPtw_Shp_Shp_Jac) + deallocate(baseline_QPtw_ShpDer_ShpDer_Jac) + deallocate(baseline_QPtw_Shp_ShpDer) + deallocate(baseline_QPtw_Shp_Jac) + deallocate(baseline_QPtw_ShpDer) + + end subroutine - ! Test and assemble variable dN*wt. - do i = 1, p%nodes_per_elem - do idx_qp = 1, p%nqp - baseline_QPtw_ShpDer(idx_qp,i) = baseline_ShpDer(i,idx_qp)*inp_QPtWeight(idx_qp) - @assertEqual(baseline_QPtw_ShpDer(idx_qp,i), p%QPtw_ShpDer(idx_qp,i), tolerance, testname) - end do - end do + @test + subroutine test_BD_InitShpDerJaco_5node() + ! branches to test + ! - 5 node, 1 element; undeformed + + ! initialize NWTC_Num constants + call SetConstants() + + tolerance = 1e-14 + + ! -------------------------------------------------------------------------- + testname = "5 node, 1 element, curved:" + + p = simpleparametertype() + p%elem_total = 1 + p%nodes_per_elem = 5 + p%nqp = 5 ! Let's use Gauss_Legendre Quadrature, which should be exact for intended polynomial test case + + ! Allocate memory for baseline results + call AllocAry(baseline_Shp , p%nodes_per_elem, p%nqp, 'Reference Shp' , ErrStat, ErrMsg) + call AllocAry(baseline_ShpDer , p%nodes_per_elem, p%nqp, 'Reference ShpDer' , ErrStat, ErrMsg) + call AllocAry(baseline_Jacobian , p%nqp, p%elem_total, 'Reference Jacobian', ErrStat, ErrMsg) + call AllocAry(baseline_QPtN , p%nqp, 'Reference QPtN' , ErrStat, ErrMsg) + call AllocAry(baseline_QPtWeight, p%nqp, 'Reference QPtWeight', ErrStat, ErrMsg) + + ! Allocate memory for other relevant variables belonging to module p + call AllocAry(baseline_QPtw_Shp_Shp_Jac , p%nqp, p%nodes_per_elem, p%nodes_per_elem, p%elem_total, 'reference QPtw_Shp_Shp_Jac' , ErrStat, ErrMsg) + call AllocAry(baseline_QPtw_ShpDer_ShpDer_Jac, p%nqp, p%nodes_per_elem, p%nodes_per_elem, p%elem_total, 'reference baseline_QPtw_ShpDer_ShpDer_Jac', ErrStat, ErrMsg) + call AllocAry(baseline_QPtw_Shp_ShpDer , p%nqp, p%nodes_per_elem, p%nodes_per_elem , 'reference QPtw_Shp_ShpDer' , ErrStat, ErrMsg) + call AllocAry(baseline_QPtw_Shp_Jac , p%nqp, p%nodes_per_elem, p%elem_total , 'reference QPtw_Shp_Jac' , ErrStat, ErrMsg) + call AllocAry(baseline_QPtw_ShpDer , p%nqp, p%nodes_per_elem , 'reference QPtw_ShpDer' , ErrStat, ErrMsg) + + ! assign baseline results + ! baseline quadrature points and weights - ! dealocate baseline variables - deallocate(baseline_Shp) - deallocate(baseline_ShpDer) - deallocate(baseline_Jacobian) - deallocate(inp_QPtWeight) - deallocate(baseline_QPtw_Shp_Shp_Jac) - deallocate(baseline_QPtw_ShpDer_ShpDer_Jac) - deallocate(baseline_QPtw_Shp_ShpDer) - deallocate(baseline_QPtw_Shp_Jac) - deallocate(baseline_QPtw_ShpDer) + baseline_QPtN(1:p%nqp) = (/ -0.9061798459386640, -0.5384693101056831, 0. , 0.5384693101056831, 0.9061798459386640 /) + baseline_QPtWeight(1:p%nqp) = (/ 0.2369268850561891, 0.4786286704993665, 0.5688888888888889, 0.4786286704993665, 0.2369268850561891 /) -end subroutine \ No newline at end of file + ! assign baseline jacobian based on example as described above + baseline_jacobian(1:p%nqp,p%elem_total) = (/ 0.6715870058501458, 1.509599209717604, 2.861380785564901, 4.097191592895223, 4.880926263217582 /) + + ! assign baseline shape functions based on example as described above + baseline_Shp(1,1:p%nqp) = (/ 0.5933706960199465, -0.10048256880508302, 0., 0.030144110771879763, -0.029205077492916114 /) + baseline_Shp(2,1:p%nqp) = (/ 0.516435198649618, 0.9313661019373962, 0., -0.09069490469997694, 0.08322282221996001 /) + baseline_Shp(3,1:p%nqp) = (/ -0.16382363939660807, 0.22966726079578503, 1., 0.22966726079578503, -0.16382363939660807 /) + baseline_Shp(4,1:p%nqp) = (/ 0.08322282221996001, -0.09069490469997694, 0., 0.9313661019373962, 0.516435198649618 /) + baseline_Shp(5,1:p%nqp) = (/ -0.029205077492916114, 0.030144110771879763, 0., -0.10048256880508302, 0.5933706960199465 /) + + ! assign baseline shape function derivatives based on example as described above + baseline_ShpDer(1,1:p%nqp) = (/ -3.705336453591454, -0.5287152679802739, 0.375, -0.24351802112960028, 0.14423640936799356 /) + baseline_ShpDer(2,1:p%nqp) = (/ 4.33282116876393, -1.0976579678283382, -1.3365845776954537, 0.7497385700132875, -0.42067623042767965 /) + baseline_ShpDer(3,1:p%nqp) = (/ -0.9039245362321631, 2.1325937846922898, 0., -2.1325937846922898, 0.9039245362321631 /) + baseline_ShpDer(4,1:p%nqp) = (/ 0.42067623042767965, -0.7497385700132875, 1.3365845776954537, 1.0976579678283382, -4.33282116876393 /) + baseline_ShpDer(5,1:p%nqp) = (/ -0.14423640936799356, 0.24351802112960028, -0.375, 0.5287152679802739, 3.705336453591454 /) + + ! Allocate memory for relevant variables belonging to module p + call AllocAry(p%Shp , p%nodes_per_elem, p%nqp, 'Shp' , ErrStat, ErrMsg) + call AllocAry(p%ShpDer , p%nodes_per_elem, p%nqp, 'ShpDer' , ErrStat, ErrMsg) + call AllocAry(p%uuN0, 3, p%nodes_per_elem, p%elem_total, 'uuN0' , ErrStat, ErrMsg) + call AllocAry(p%Jacobian , p%elem_total , p%nqp, 'Jacobian' , ErrStat, ErrMsg) ! reversed arguments? + call AllocAry(p%QPtN , p%nqp , 'QPtN' , ErrStat, ErrMsg) + call AllocAry(p%QPtWeight, p%nqp , 'QPtWeight', ErrStat, ErrMsg) + + ! Allocate memory for other relevant variables belonging to module p + call AllocAry(p%QPtw_Shp_Shp_Jac , p%nqp, p%nodes_per_elem, p%nodes_per_elem,p%elem_total, 'QPtw_Shp_Shp_Jac', ErrStat, ErrMsg) + call AllocAry(p%QPtw_Shp_ShpDer , p%nqp, p%nodes_per_elem, p%nodes_per_elem , 'QPtw_Shp_ShpDer', ErrStat, ErrMsg) + call AllocAry(p%QPtw_ShpDer_ShpDer_Jac, p%nqp, p%nodes_per_elem, p%nodes_per_elem,p%elem_total, 'QPtw_ShpDer_ShpDer_Jac', ErrStat, ErrMsg) + call AllocAry(p%QPtw_Shp_Jac , p%nqp, p%nodes_per_elem, p%elem_total , 'QPtw_Shp_Jac', ErrStat, ErrMsg) + call AllocAry(p%QPtw_ShpDer , p%nqp, p%nodes_per_elem , 'QPtw_ShpDer', ErrStat, ErrMsg) + + ! uuN0 is of dimension (3 dof, nodes_per_elem, elem_total) + p%uuN0(1:3,1,1) = (/ 0.0, 0.0, 0.0 /) + p%uuN0(1:3,2,1) = (/ 0.16237631096713473, 0.17578464768961147, 0.1481911137890286 /) + p%uuN0(1:3,3,1) = (/ 0.25, 1., 1.1875 /) + p%uuN0(1:3,4,1) = (/ -0.30523345382427747, 2.4670724951675314, 2.953849702537502 /) + p%uuN0(1:3,5,1) = (/ -1., 3.5, 4. /) + + !p%QPtN = (/ 0.5 /) ! Note, we assume 1 quadrature point + !p%QPtWeight = inp_QPtWeight ! Since we assume 1 quadrature point, the weight by defualt = 1 + + ! Using BD_GaussPointWeight; hoping it's tested! + call BD_GaussPointWeight(p%nqp, p%QPtN, p%QPtWeight, ErrStat, ErrMsg) + + @assertEqual(baseline_QPtN, p%QPtN , tolerance, testname) + @assertEqual(baseline_QPtWeight, p%QPtWeight, tolerance, testname) + + ! Allocate memory for GLL node positions in 1D parametric space +! call AllocAry(gll_nodes, p%nodes_per_elem, "GLL points array", ErrStat, ErrMsg) +! gll_nodes = (/ -1.,-0.6546536707079771,0.,0.6546536707079771,1. /) +! +! ! call the test subroutine +! call BD_InitShpDerJaco(gll_nodes, p) +! +! ! check the baseline shape functions and their derivatives +! do idx_qp = 1, p%nqp +! do j = 1, p%nodes_per_elem +! @assertEqual(baseline_Shp(j,idx_qp) , p%Shp(j,idx_qp) , tolerance, testname) +! @assertEqual(baseline_ShpDer(j,idx_qp), p%ShpDer(j,idx_qp), tolerance, testname) +! end do +! end do +! +! ! check the baseline jacobian +! do nelem = 1, p%elem_total +! do idx_qp = 1, p%nqp +! @assertEqual(baseline_jacobian(nelem,idx_qp), p%jacobian(nelem,idx_qp), tolerance, testname) +! end do +! end do +! +! ! Test and assemble variables N*N^T*wt*Jacobian and dN*dN^T*wt/Jacobian +! do nelem = 1, p%elem_total +! do idx_qp = 1, p%nqp +! do j = 1, p%nodes_per_elem +! do i = 1, p%nodes_per_elem +! ! Check the variable N*N^T*Jacobian +! baseline_QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem) = baseline_Shp(i,idx_qp)*baseline_Shp(j,idx_qp)*baseline_QPtWeight(idx_qp)*baseline_jacobian(idx_qp,nelem) +! @assertEqual(baseline_QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem), p%QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem), tolerance, testname) +! +! ! Check the variable dN*dN^T*Jacobian +! baseline_QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem) = baseline_ShpDer(i,idx_qp)*baseline_ShpDer(j,idx_qp)*baseline_QPtWeight(idx_qp)/baseline_jacobian(idx_qp,nelem) +! @assertEqual(baseline_QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem), p%QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem), tolerance, testname) +! end do +! end do +! end do +! end do +! +! ! Test and assemble variable N*dN^T*wt*Jacobian +! do idx_qp = 1, p%nqp +! do j = 1, p%nodes_per_elem +! do i = 1, p%nodes_per_elem +! baseline_QPtw_Shp_ShpDer(idx_qp,i,j) = baseline_Shp(i,idx_qp)*baseline_ShpDer(j,idx_qp)*baseline_QPtWeight(idx_qp) +! @assertEqual(baseline_QPtw_Shp_ShpDer(idx_qp,i,j), p%QPtw_Shp_ShpDer(idx_qp,i,j), tolerance, testname) +! end do +! end do +! end do +! +! ! Test and assemble variable N*wt*Jacobian +! do nelem = 1, p%elem_total +! do i = 1, p%nodes_per_elem +! do idx_qp = 1, p%nqp +! baseline_QPtw_Shp_Jac(idx_qp,i,nelem) = baseline_Shp(i,idx_qp)*baseline_QPtWeight(idx_qp)*baseline_Jacobian(idx_qp,nelem) +! @assertEqual(baseline_QPtw_Shp_Jac(idx_qp,i,nelem), p%QPtw_Shp_Jac(idx_qp,i,nelem), tolerance, testname) +! end do +! end do +! end do +! +! ! Test and assemble variable dN*wt. +! do i = 1, p%nodes_per_elem +! do idx_qp = 1, p%nqp +! baseline_QPtw_ShpDer(idx_qp,i) = baseline_ShpDer(i,idx_qp)*baseline_QPtWeight(idx_qp) +! @assertEqual(baseline_QPtw_ShpDer(idx_qp,i), p%QPtw_ShpDer(idx_qp,i), tolerance, testname) +! end do +! end do +! + ! dealocate baseline variables + deallocate(baseline_Shp) + deallocate(baseline_ShpDer) + deallocate(baseline_Jacobian) + deallocate(baseline_QPtN) + deallocate(baseline_QPtWeight) + deallocate(baseline_QPtw_Shp_Shp_Jac) + deallocate(baseline_QPtw_ShpDer_ShpDer_Jac) + deallocate(baseline_QPtw_Shp_ShpDer) + deallocate(baseline_QPtw_Shp_Jac) + deallocate(baseline_QPtw_ShpDer) + + end subroutine +end module From 982b24652db773046d04ac153b9592db452cc1d8 Mon Sep 17 00:00:00 2001 From: Mike Sprague Date: Thu, 8 Oct 2020 15:22:45 -0600 Subject: [PATCH 3/8] added simpleparametertype_teardown; reworkd test_BD_InitShpDerJaco.F90 for a 5-node curved element --- .../beamdyn/tests/test_BD_DistrLoadCopy.F90 | 4 +- .../beamdyn/tests/test_BD_GravityForce.F90 | 6 +- .../beamdyn/tests/test_BD_InitShpDerJaco.F90 | 292 +++--------------- .../beamdyn/tests/test_BD_QPData_mEta_rho.F90 | 3 +- modules/beamdyn/tests/test_BD_diffmtc.F90 | 18 +- .../tests/test_ExtractRelativeRotation.F90 | 2 +- modules/beamdyn/tests/test_tools.F90 | 48 ++- 7 files changed, 111 insertions(+), 262 deletions(-) diff --git a/modules/beamdyn/tests/test_BD_DistrLoadCopy.F90 b/modules/beamdyn/tests/test_BD_DistrLoadCopy.F90 index 1e21931dd3..88e5ae87f6 100644 --- a/modules/beamdyn/tests/test_BD_DistrLoadCopy.F90 +++ b/modules/beamdyn/tests/test_BD_DistrLoadCopy.F90 @@ -29,7 +29,7 @@ subroutine test_BD_DistrLoadCopy() testname = "static simple beam under gravity:" ! build the parametertype, inputtype, and miscvartype - parametertype = simpleParameterType() + parametertype = simpleParameterType(1,16,16,0) miscvartype = simpleMiscVarType(parametertype%nqp, parametertype%elem_total) inputtype = simpleInputType(parametertype%nqp, parametertype%elem_total) @@ -41,4 +41,6 @@ subroutine test_BD_DistrLoadCopy() @assertEqual((/ -9*(j-1)-3*(i-1)-1, -9*(j-1)-3*(i-1)-2, -9*(j-1)-3*(i-1)-3 /), miscvartype%DistrLoad_QP(4:6,i,j)) end do end do + + call simpleparametertype_teardown(parametertype) end subroutine diff --git a/modules/beamdyn/tests/test_BD_GravityForce.F90 b/modules/beamdyn/tests/test_BD_GravityForce.F90 index 3e94d21368..aa99473122 100644 --- a/modules/beamdyn/tests/test_BD_GravityForce.F90 +++ b/modules/beamdyn/tests/test_BD_GravityForce.F90 @@ -30,7 +30,7 @@ subroutine test_BD_GravityForce() baseline(4:6) = (/ 0.0, 0.0, 0.0 /) ! allocate and build the custom types - parametertype = simpleParameterType() + parametertype = simpleParameterType(1,16,16,0) miscvartype = simpleMiscVarType(parametertype%nqp, parametertype%elem_total) gravity = getGravityInZ() @@ -40,5 +40,7 @@ subroutine test_BD_GravityForce() ! test the values @assertEqual(baseline, miscvartype%qp%Fg(:,1,1), tolerance, testname) - + + call simpleparametertype_teardown(parametertype) + end subroutine diff --git a/modules/beamdyn/tests/test_BD_InitShpDerJaco.F90 b/modules/beamdyn/tests/test_BD_InitShpDerJaco.F90 index cfba996cbf..a1e88cb038 100644 --- a/modules/beamdyn/tests/test_BD_InitShpDerJaco.F90 +++ b/modules/beamdyn/tests/test_BD_InitShpDerJaco.F90 @@ -22,51 +22,29 @@ module test_BD_InitShpDerJaco contains @test - subroutine test_BD_InitShpDerJaco_3node() + subroutine test_BD_InitShpDerJaco_5node() ! branches to test - ! - 3 node, 1 element; undeformed - + ! - 5 node, 1 element; undeformed ! initialize NWTC_Num constants call SetConstants() - tolerance = 1e-14 - + tolerance = 1e-13 ! -------------------------------------------------------------------------- - testname = "3 node, 1 element, undeformed:" - - ! Lets assume a 3 node element with parametric coordinate s where s lies in [-1 1] - ! The three nodes lie at s1 = -1, s2 = 0, and s3 = 1 - ! The corresponding positions of the nodes in physical space are assumed to be - ! x1 = 0, x2 = 10, and x2 = 20 - ! - ! The quadratic shape functions for such a 3 node element at any given 's' are given by - ! N1(s) = (s2-s)*(s3-s)/(s2-s1)/(s3-s1) - ! N2(s) = (s1-s)*(s3-s)/(s1-s2)/(s3-s2) - ! N3(s) = (s1-s)*(s2-s)/(s1-s3)/(s2-s3) - ! - ! We are interested in testing the function test_BD_InitShpDerJaco for s = 0.5 - ! At s = 0.5, the above values are computed to be N1 = -0.125, N2 = 0.75, and N3 = -0.375 - ! - ! The derivatives dN/ds are computed to be dN1 = -, dN2 = -1, and dN3 = 1 - ! - ! the Jacobian is invariant of the quadrature point and account for change in volume - ! between physical and parametric space. Jacobian = meas(x)/meas(s), where meas is the - ! measured quantity (volume in 3D, area in 2D, and lenght in 1D) - ! The Jacobian for this element is given by 20/2 = 10 + testname = "5 node, 1 element, curved:" - ! build the p object based on the above mentioned test model - p = simpleparametertype() - p%elem_total = 1 - p%nodes_per_elem = 3 - p%nqp = 1 + p = simpleparametertype(1,5,5,0) + !p%elem_total = 1 + !p%nodes_per_elem = 5 + !p%nqp = 5 ! Let's use Gauss_Legendre Quadrature, which should be exact for intended polynomial test case - ! Allocate memory for baseline results FIXME -- seems like arguments reversed + ! Allocate memory for baseline results call AllocAry(baseline_Shp , p%nodes_per_elem, p%nqp, 'Reference Shp' , ErrStat, ErrMsg) call AllocAry(baseline_ShpDer , p%nodes_per_elem, p%nqp, 'Reference ShpDer' , ErrStat, ErrMsg) - call AllocAry(baseline_Jacobian, p%elem_total , p%nqp, 'Reference Jacobian', ErrStat, ErrMsg) - call AllocAry(inp_QPtWeight , p%nqp , 'QPtWeight' , ErrStat, ErrMsg) + call AllocAry(baseline_Jacobian , p%nqp, p%elem_total, 'Reference Jacobian', ErrStat, ErrMsg) + call AllocAry(baseline_QPtN , p%nqp, 'Reference QPtN' , ErrStat, ErrMsg) + call AllocAry(baseline_QPtWeight, p%nqp, 'Reference QPtWeight', ErrStat, ErrMsg) ! Allocate memory for other relevant variables belonging to module p call AllocAry(baseline_QPtw_Shp_Shp_Jac , p%nqp, p%nodes_per_elem, p%nodes_per_elem, p%elem_total, 'reference QPtw_Shp_Shp_Jac' , ErrStat, ErrMsg) @@ -75,50 +53,45 @@ subroutine test_BD_InitShpDerJaco_3node() call AllocAry(baseline_QPtw_Shp_Jac , p%nqp, p%nodes_per_elem, p%elem_total , 'reference QPtw_Shp_Jac' , ErrStat, ErrMsg) call AllocAry(baseline_QPtw_ShpDer , p%nqp, p%nodes_per_elem , 'reference QPtw_ShpDer' , ErrStat, ErrMsg) - ! assign baseling results - ! assign baseline jacobian based on example as described above - baseline_jacobian(1,1) = 10.0 ! we assume 1 element quadrature point. Hence we set only the index (1,1) + ! assign baseline results + ! baseline quadrature points and weights; this is 5-point Gauss-Legendre quadrature + + baseline_QPtN(1:p%nqp) = (/ -0.9061798459386640, -0.5384693101056831, 0. , 0.5384693101056831, 0.9061798459386640 /) + baseline_QPtWeight(1:p%nqp) = (/ 0.2369268850561891, 0.4786286704993665, 0.5688888888888889, 0.4786286704993665, 0.2369268850561891 /) + + ! assign baseline jacobian based; these values were calculated in separte mathematica script + baseline_jacobian(1:p%nqp,1) = (/ 0.6715870058501458, 1.509599209717604, 2.861380785564901, 4.097191592895223, 4.880926263217582 /) ! assign baseline shape functions based on example as described above - baseline_Shp(1,1) = -0.125 - baseline_Shp(2,1) = 0.750 - baseline_Shp(3,1) = 0.375 + baseline_Shp(1,1:p%nqp) = (/ 0.5933706960199465, -0.10048256880508302, 0., 0.030144110771879763, -0.029205077492916114 /) + baseline_Shp(2,1:p%nqp) = (/ 0.516435198649618, 0.9313661019373962, 0., -0.09069490469997694, 0.08322282221996001 /) + baseline_Shp(3,1:p%nqp) = (/ -0.16382363939660807, 0.22966726079578503, 1., 0.22966726079578503, -0.16382363939660807 /) + baseline_Shp(4,1:p%nqp) = (/ 0.08322282221996001, -0.09069490469997694, 0., 0.9313661019373962, 0.516435198649618 /) + baseline_Shp(5,1:p%nqp) = (/ -0.029205077492916114, 0.030144110771879763, 0., -0.10048256880508302, 0.5933706960199465 /) ! assign baseline shape function derivatives based on example as described above - baseline_ShpDer(1,1) = 0.0 - baseline_ShpDer(2,1) = -1.0 - baseline_ShpDer(3,1) = 1.0 - - ! assign the weight to the quadrature point which is an input parameter - inp_QPtWeight(1) = 1.0; - - ! Allocate memory for relevant variables belonging to module p FIXME -- seems like arguments reversed; work becuase nodes_per_elem = nqp = 1 - call AllocAry(p%Shp , p%nodes_per_elem, p%nqp, 'Shp' , ErrStat, ErrMsg) - call AllocAry(p%ShpDer , p%nodes_per_elem, p%nqp, 'ShpDer' , ErrStat, ErrMsg) - call AllocAry(p%uuN0, 3, p%nodes_per_elem, p%elem_total, 'uuN0' , ErrStat, ErrMsg) - call AllocAry(p%Jacobian , p%elem_total , p%nqp, 'Jacobian' , ErrStat, ErrMsg) ! reversed arguments? - call AllocAry(p%QPtN , p%nodes_per_elem , 'QPtN' , ErrStat, ErrMsg) - call AllocAry(p%QPtWeight, p%nqp , 'QPtWeight', ErrStat, ErrMsg) - - ! Allocate memory for other relevant variables belonging to module p - call AllocAry(p%QPtw_Shp_Shp_Jac , p%nqp, p%nodes_per_elem, p%nodes_per_elem,p%elem_total, 'QPtw_Shp_Shp_Jac', ErrStat, ErrMsg) - call AllocAry(p%QPtw_Shp_ShpDer , p%nqp, p%nodes_per_elem, p%nodes_per_elem , 'QPtw_Shp_ShpDer', ErrStat, ErrMsg) - call AllocAry(p%QPtw_ShpDer_ShpDer_Jac, p%nqp, p%nodes_per_elem, p%nodes_per_elem,p%elem_total, 'QPtw_ShpDer_ShpDer_Jac', ErrStat, ErrMsg) - call AllocAry(p%QPtw_Shp_Jac , p%nqp, p%nodes_per_elem, p%elem_total , 'QPtw_Shp_Jac', ErrStat, ErrMsg) - call AllocAry(p%QPtw_ShpDer , p%nqp, p%nodes_per_elem , 'QPtw_ShpDer', ErrStat, ErrMsg) + baseline_ShpDer(1,1:p%nqp) = (/ -3.705336453591454, -0.5287152679802739, 0.375, -0.24351802112960028, 0.14423640936799356 /) + baseline_ShpDer(2,1:p%nqp) = (/ 4.33282116876393, -1.0976579678283382, -1.3365845776954537, 0.7497385700132875, -0.42067623042767965 /) + baseline_ShpDer(3,1:p%nqp) = (/ -0.9039245362321631, 2.1325937846922898, 0., -2.1325937846922898, 0.9039245362321631 /) + baseline_ShpDer(4,1:p%nqp) = (/ 0.42067623042767965, -0.7497385700132875, 1.3365845776954537, 1.0976579678283382, -4.33282116876393 /) + baseline_ShpDer(5,1:p%nqp) = (/ -0.14423640936799356, 0.24351802112960028, -0.375, 0.5287152679802739, 3.705336453591454 /) ! uuN0 is of dimension (3 dof, nodes_per_elem, elem_total) p%uuN0(1:3,1,1) = (/ 0.0, 0.0, 0.0 /) - p%uuN0(1:3,2,1) = (/ 0.0, 0.0, 10.0 /) - p%uuN0(1:3,3,1) = (/ 0.0, 0.0, 20.0 /) - - p%QPtN = (/ 0.5 /) ! Note, we assume 1 quadrature point - - p%QPtWeight = inp_QPtWeight ! Since we assume 1 quadrature point, the weight by defualt = 1 - + p%uuN0(1:3,2,1) = (/ 0.16237631096713473, 0.17578464768961147, 0.1481911137890286 /) + p%uuN0(1:3,3,1) = (/ 0.25, 1., 1.1875 /) + p%uuN0(1:3,4,1) = (/ -0.30523345382427747, 2.4670724951675314, 2.953849702537502 /) + p%uuN0(1:3,5,1) = (/ -1., 3.5, 4. /) + + ! Using BD_GaussPointWeight; hoping it's tested! + call BD_GaussPointWeight(p%nqp, p%QPtN, p%QPtWeight, ErrStat, ErrMsg) + + @assertEqual(baseline_QPtN, p%QPtN , tolerance, testname) + @assertEqual(baseline_QPtWeight, p%QPtWeight, tolerance, testname) + ! Allocate memory for GLL node positions in 1D parametric space call AllocAry(gll_nodes, p%nodes_per_elem, "GLL points array", ErrStat, ErrMsg) - gll_nodes = (/ -1.0, 0.0, 1.0 /) + gll_nodes = (/ -1., -0.6546536707079771, 0., 0.6546536707079771, 1. /) ! call the test subroutine call BD_InitShpDerJaco(gll_nodes, p) @@ -134,7 +107,7 @@ subroutine test_BD_InitShpDerJaco_3node() ! check the baseline jacobian do nelem = 1, p%elem_total do idx_qp = 1, p%nqp - @assertEqual(baseline_jacobian(nelem,idx_qp), p%jacobian(nelem,idx_qp), tolerance, testname) + @assertEqual(baseline_jacobian(idx_qp,nelem), p%jacobian(idx_qp,nelem), tolerance, testname) end do end do @@ -144,11 +117,11 @@ subroutine test_BD_InitShpDerJaco_3node() do j = 1, p%nodes_per_elem do i = 1, p%nodes_per_elem ! Check the variable N*N^T*Jacobian - baseline_QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem) = baseline_Shp(i,idx_qp)*baseline_Shp(j,idx_qp)*inp_QPtWeight(idx_qp)*baseline_jacobian(idx_qp,nelem) + baseline_QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem) = baseline_Shp(i,idx_qp)*baseline_Shp(j,idx_qp)*baseline_QPtWeight(idx_qp)*baseline_jacobian(idx_qp,nelem) @assertEqual(baseline_QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem), p%QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem), tolerance, testname) ! Check the variable dN*dN^T*Jacobian - baseline_QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem) = baseline_ShpDer(i,idx_qp)*baseline_ShpDer(j,idx_qp)*inp_QPtWeight(idx_qp)/baseline_jacobian(idx_qp,nelem) + baseline_QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem) = baseline_ShpDer(i,idx_qp)*baseline_ShpDer(j,idx_qp)*baseline_QPtWeight(idx_qp)/baseline_jacobian(idx_qp,nelem) @assertEqual(baseline_QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem), p%QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem), tolerance, testname) end do end do @@ -159,7 +132,7 @@ subroutine test_BD_InitShpDerJaco_3node() do idx_qp = 1, p%nqp do j = 1, p%nodes_per_elem do i = 1, p%nodes_per_elem - baseline_QPtw_Shp_ShpDer(idx_qp,i,j) = baseline_Shp(i,idx_qp)*baseline_ShpDer(j,idx_qp)*inp_QPtWeight(idx_qp) + baseline_QPtw_Shp_ShpDer(idx_qp,i,j) = baseline_Shp(i,idx_qp)*baseline_ShpDer(j,idx_qp)*baseline_QPtWeight(idx_qp) @assertEqual(baseline_QPtw_Shp_ShpDer(idx_qp,i,j), p%QPtw_Shp_ShpDer(idx_qp,i,j), tolerance, testname) end do end do @@ -169,7 +142,7 @@ subroutine test_BD_InitShpDerJaco_3node() do nelem = 1, p%elem_total do i = 1, p%nodes_per_elem do idx_qp = 1, p%nqp - baseline_QPtw_Shp_Jac(idx_qp,i,nelem) = baseline_Shp(i,idx_qp)*inp_QPtWeight(idx_qp)*baseline_Jacobian(idx_qp,nelem) + baseline_QPtw_Shp_Jac(idx_qp,i,nelem) = baseline_Shp(i,idx_qp)*baseline_QPtWeight(idx_qp)*baseline_Jacobian(idx_qp,nelem) @assertEqual(baseline_QPtw_Shp_Jac(idx_qp,i,nelem), p%QPtw_Shp_Jac(idx_qp,i,nelem), tolerance, testname) end do end do @@ -178,178 +151,13 @@ subroutine test_BD_InitShpDerJaco_3node() ! Test and assemble variable dN*wt. do i = 1, p%nodes_per_elem do idx_qp = 1, p%nqp - baseline_QPtw_ShpDer(idx_qp,i) = baseline_ShpDer(i,idx_qp)*inp_QPtWeight(idx_qp) + baseline_QPtw_ShpDer(idx_qp,i) = baseline_ShpDer(i,idx_qp)*baseline_QPtWeight(idx_qp) @assertEqual(baseline_QPtw_ShpDer(idx_qp,i), p%QPtw_ShpDer(idx_qp,i), tolerance, testname) end do end do ! dealocate baseline variables - deallocate(baseline_Shp) - deallocate(baseline_ShpDer) - deallocate(baseline_Jacobian) - deallocate(inp_QPtWeight) - deallocate(baseline_QPtw_Shp_Shp_Jac) - deallocate(baseline_QPtw_ShpDer_ShpDer_Jac) - deallocate(baseline_QPtw_Shp_ShpDer) - deallocate(baseline_QPtw_Shp_Jac) - deallocate(baseline_QPtw_ShpDer) - - end subroutine - - @test - subroutine test_BD_InitShpDerJaco_5node() - ! branches to test - ! - 5 node, 1 element; undeformed - - ! initialize NWTC_Num constants - call SetConstants() - - tolerance = 1e-14 - - ! -------------------------------------------------------------------------- - testname = "5 node, 1 element, curved:" - - p = simpleparametertype() - p%elem_total = 1 - p%nodes_per_elem = 5 - p%nqp = 5 ! Let's use Gauss_Legendre Quadrature, which should be exact for intended polynomial test case - - ! Allocate memory for baseline results - call AllocAry(baseline_Shp , p%nodes_per_elem, p%nqp, 'Reference Shp' , ErrStat, ErrMsg) - call AllocAry(baseline_ShpDer , p%nodes_per_elem, p%nqp, 'Reference ShpDer' , ErrStat, ErrMsg) - call AllocAry(baseline_Jacobian , p%nqp, p%elem_total, 'Reference Jacobian', ErrStat, ErrMsg) - call AllocAry(baseline_QPtN , p%nqp, 'Reference QPtN' , ErrStat, ErrMsg) - call AllocAry(baseline_QPtWeight, p%nqp, 'Reference QPtWeight', ErrStat, ErrMsg) - - ! Allocate memory for other relevant variables belonging to module p - call AllocAry(baseline_QPtw_Shp_Shp_Jac , p%nqp, p%nodes_per_elem, p%nodes_per_elem, p%elem_total, 'reference QPtw_Shp_Shp_Jac' , ErrStat, ErrMsg) - call AllocAry(baseline_QPtw_ShpDer_ShpDer_Jac, p%nqp, p%nodes_per_elem, p%nodes_per_elem, p%elem_total, 'reference baseline_QPtw_ShpDer_ShpDer_Jac', ErrStat, ErrMsg) - call AllocAry(baseline_QPtw_Shp_ShpDer , p%nqp, p%nodes_per_elem, p%nodes_per_elem , 'reference QPtw_Shp_ShpDer' , ErrStat, ErrMsg) - call AllocAry(baseline_QPtw_Shp_Jac , p%nqp, p%nodes_per_elem, p%elem_total , 'reference QPtw_Shp_Jac' , ErrStat, ErrMsg) - call AllocAry(baseline_QPtw_ShpDer , p%nqp, p%nodes_per_elem , 'reference QPtw_ShpDer' , ErrStat, ErrMsg) - - ! assign baseline results - ! baseline quadrature points and weights - - baseline_QPtN(1:p%nqp) = (/ -0.9061798459386640, -0.5384693101056831, 0. , 0.5384693101056831, 0.9061798459386640 /) - baseline_QPtWeight(1:p%nqp) = (/ 0.2369268850561891, 0.4786286704993665, 0.5688888888888889, 0.4786286704993665, 0.2369268850561891 /) - - ! assign baseline jacobian based on example as described above - baseline_jacobian(1:p%nqp,p%elem_total) = (/ 0.6715870058501458, 1.509599209717604, 2.861380785564901, 4.097191592895223, 4.880926263217582 /) - - ! assign baseline shape functions based on example as described above - baseline_Shp(1,1:p%nqp) = (/ 0.5933706960199465, -0.10048256880508302, 0., 0.030144110771879763, -0.029205077492916114 /) - baseline_Shp(2,1:p%nqp) = (/ 0.516435198649618, 0.9313661019373962, 0., -0.09069490469997694, 0.08322282221996001 /) - baseline_Shp(3,1:p%nqp) = (/ -0.16382363939660807, 0.22966726079578503, 1., 0.22966726079578503, -0.16382363939660807 /) - baseline_Shp(4,1:p%nqp) = (/ 0.08322282221996001, -0.09069490469997694, 0., 0.9313661019373962, 0.516435198649618 /) - baseline_Shp(5,1:p%nqp) = (/ -0.029205077492916114, 0.030144110771879763, 0., -0.10048256880508302, 0.5933706960199465 /) - - ! assign baseline shape function derivatives based on example as described above - baseline_ShpDer(1,1:p%nqp) = (/ -3.705336453591454, -0.5287152679802739, 0.375, -0.24351802112960028, 0.14423640936799356 /) - baseline_ShpDer(2,1:p%nqp) = (/ 4.33282116876393, -1.0976579678283382, -1.3365845776954537, 0.7497385700132875, -0.42067623042767965 /) - baseline_ShpDer(3,1:p%nqp) = (/ -0.9039245362321631, 2.1325937846922898, 0., -2.1325937846922898, 0.9039245362321631 /) - baseline_ShpDer(4,1:p%nqp) = (/ 0.42067623042767965, -0.7497385700132875, 1.3365845776954537, 1.0976579678283382, -4.33282116876393 /) - baseline_ShpDer(5,1:p%nqp) = (/ -0.14423640936799356, 0.24351802112960028, -0.375, 0.5287152679802739, 3.705336453591454 /) - - ! Allocate memory for relevant variables belonging to module p - call AllocAry(p%Shp , p%nodes_per_elem, p%nqp, 'Shp' , ErrStat, ErrMsg) - call AllocAry(p%ShpDer , p%nodes_per_elem, p%nqp, 'ShpDer' , ErrStat, ErrMsg) - call AllocAry(p%uuN0, 3, p%nodes_per_elem, p%elem_total, 'uuN0' , ErrStat, ErrMsg) - call AllocAry(p%Jacobian , p%elem_total , p%nqp, 'Jacobian' , ErrStat, ErrMsg) ! reversed arguments? - call AllocAry(p%QPtN , p%nqp , 'QPtN' , ErrStat, ErrMsg) - call AllocAry(p%QPtWeight, p%nqp , 'QPtWeight', ErrStat, ErrMsg) - - ! Allocate memory for other relevant variables belonging to module p - call AllocAry(p%QPtw_Shp_Shp_Jac , p%nqp, p%nodes_per_elem, p%nodes_per_elem,p%elem_total, 'QPtw_Shp_Shp_Jac', ErrStat, ErrMsg) - call AllocAry(p%QPtw_Shp_ShpDer , p%nqp, p%nodes_per_elem, p%nodes_per_elem , 'QPtw_Shp_ShpDer', ErrStat, ErrMsg) - call AllocAry(p%QPtw_ShpDer_ShpDer_Jac, p%nqp, p%nodes_per_elem, p%nodes_per_elem,p%elem_total, 'QPtw_ShpDer_ShpDer_Jac', ErrStat, ErrMsg) - call AllocAry(p%QPtw_Shp_Jac , p%nqp, p%nodes_per_elem, p%elem_total , 'QPtw_Shp_Jac', ErrStat, ErrMsg) - call AllocAry(p%QPtw_ShpDer , p%nqp, p%nodes_per_elem , 'QPtw_ShpDer', ErrStat, ErrMsg) - - ! uuN0 is of dimension (3 dof, nodes_per_elem, elem_total) - p%uuN0(1:3,1,1) = (/ 0.0, 0.0, 0.0 /) - p%uuN0(1:3,2,1) = (/ 0.16237631096713473, 0.17578464768961147, 0.1481911137890286 /) - p%uuN0(1:3,3,1) = (/ 0.25, 1., 1.1875 /) - p%uuN0(1:3,4,1) = (/ -0.30523345382427747, 2.4670724951675314, 2.953849702537502 /) - p%uuN0(1:3,5,1) = (/ -1., 3.5, 4. /) - - !p%QPtN = (/ 0.5 /) ! Note, we assume 1 quadrature point - !p%QPtWeight = inp_QPtWeight ! Since we assume 1 quadrature point, the weight by defualt = 1 - - ! Using BD_GaussPointWeight; hoping it's tested! - call BD_GaussPointWeight(p%nqp, p%QPtN, p%QPtWeight, ErrStat, ErrMsg) - - @assertEqual(baseline_QPtN, p%QPtN , tolerance, testname) - @assertEqual(baseline_QPtWeight, p%QPtWeight, tolerance, testname) - - ! Allocate memory for GLL node positions in 1D parametric space -! call AllocAry(gll_nodes, p%nodes_per_elem, "GLL points array", ErrStat, ErrMsg) -! gll_nodes = (/ -1.,-0.6546536707079771,0.,0.6546536707079771,1. /) -! -! ! call the test subroutine -! call BD_InitShpDerJaco(gll_nodes, p) -! -! ! check the baseline shape functions and their derivatives -! do idx_qp = 1, p%nqp -! do j = 1, p%nodes_per_elem -! @assertEqual(baseline_Shp(j,idx_qp) , p%Shp(j,idx_qp) , tolerance, testname) -! @assertEqual(baseline_ShpDer(j,idx_qp), p%ShpDer(j,idx_qp), tolerance, testname) -! end do -! end do -! -! ! check the baseline jacobian -! do nelem = 1, p%elem_total -! do idx_qp = 1, p%nqp -! @assertEqual(baseline_jacobian(nelem,idx_qp), p%jacobian(nelem,idx_qp), tolerance, testname) -! end do -! end do -! -! ! Test and assemble variables N*N^T*wt*Jacobian and dN*dN^T*wt/Jacobian -! do nelem = 1, p%elem_total -! do idx_qp = 1, p%nqp -! do j = 1, p%nodes_per_elem -! do i = 1, p%nodes_per_elem -! ! Check the variable N*N^T*Jacobian -! baseline_QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem) = baseline_Shp(i,idx_qp)*baseline_Shp(j,idx_qp)*baseline_QPtWeight(idx_qp)*baseline_jacobian(idx_qp,nelem) -! @assertEqual(baseline_QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem), p%QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem), tolerance, testname) -! -! ! Check the variable dN*dN^T*Jacobian -! baseline_QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem) = baseline_ShpDer(i,idx_qp)*baseline_ShpDer(j,idx_qp)*baseline_QPtWeight(idx_qp)/baseline_jacobian(idx_qp,nelem) -! @assertEqual(baseline_QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem), p%QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem), tolerance, testname) -! end do -! end do -! end do -! end do -! -! ! Test and assemble variable N*dN^T*wt*Jacobian -! do idx_qp = 1, p%nqp -! do j = 1, p%nodes_per_elem -! do i = 1, p%nodes_per_elem -! baseline_QPtw_Shp_ShpDer(idx_qp,i,j) = baseline_Shp(i,idx_qp)*baseline_ShpDer(j,idx_qp)*baseline_QPtWeight(idx_qp) -! @assertEqual(baseline_QPtw_Shp_ShpDer(idx_qp,i,j), p%QPtw_Shp_ShpDer(idx_qp,i,j), tolerance, testname) -! end do -! end do -! end do -! -! ! Test and assemble variable N*wt*Jacobian -! do nelem = 1, p%elem_total -! do i = 1, p%nodes_per_elem -! do idx_qp = 1, p%nqp -! baseline_QPtw_Shp_Jac(idx_qp,i,nelem) = baseline_Shp(i,idx_qp)*baseline_QPtWeight(idx_qp)*baseline_Jacobian(idx_qp,nelem) -! @assertEqual(baseline_QPtw_Shp_Jac(idx_qp,i,nelem), p%QPtw_Shp_Jac(idx_qp,i,nelem), tolerance, testname) -! end do -! end do -! end do -! -! ! Test and assemble variable dN*wt. -! do i = 1, p%nodes_per_elem -! do idx_qp = 1, p%nqp -! baseline_QPtw_ShpDer(idx_qp,i) = baseline_ShpDer(i,idx_qp)*baseline_QPtWeight(idx_qp) -! @assertEqual(baseline_QPtw_ShpDer(idx_qp,i), p%QPtw_ShpDer(idx_qp,i), tolerance, testname) -! end do -! end do -! - ! dealocate baseline variables + if (allocated(gll_nodes)) deallocate(gll_nodes) deallocate(baseline_Shp) deallocate(baseline_ShpDer) deallocate(baseline_Jacobian) @@ -360,6 +168,8 @@ subroutine test_BD_InitShpDerJaco_5node() deallocate(baseline_QPtw_Shp_ShpDer) deallocate(baseline_QPtw_Shp_Jac) deallocate(baseline_QPtw_ShpDer) + + call SimpleParameterType_TearDown(p) end subroutine end module diff --git a/modules/beamdyn/tests/test_BD_QPData_mEta_rho.F90 b/modules/beamdyn/tests/test_BD_QPData_mEta_rho.F90 index ca34cf1477..3382904b2d 100644 --- a/modules/beamdyn/tests/test_BD_QPData_mEta_rho.F90 +++ b/modules/beamdyn/tests/test_BD_QPData_mEta_rho.F90 @@ -34,7 +34,7 @@ subroutine test_BD_QPData_mEta_rho() baselineRR0mEta = (/ 0.0, 0.0, 0.0 /) ! allocate and build the custom input types - parametertype = simpleParameterType() + parametertype = simpleParameterType(1,16,16,0) miscvartype = simpleMiscVarType(parametertype%nqp, parametertype%elem_total) ! allocate the results @@ -46,4 +46,5 @@ subroutine test_BD_QPData_mEta_rho() @assertEqual(baselineRR0mEta, miscvartype%qp%RR0mEta(:,i,j), tolerance, testname) end do end do + call simpleParameterType_teardown(parametertype) end subroutine diff --git a/modules/beamdyn/tests/test_BD_diffmtc.F90 b/modules/beamdyn/tests/test_BD_diffmtc.F90 index ca5810d26a..84093e5423 100644 --- a/modules/beamdyn/tests/test_BD_diffmtc.F90 +++ b/modules/beamdyn/tests/test_BD_diffmtc.F90 @@ -68,9 +68,9 @@ subroutine test_BD_diffmtc_2node() ! shpder(1,:) = -0.5, -0.5 ! shpder(2,:) = 0.5, 0.5 - parametertype = simpleParameterType() - parametertype%nodes_per_elem = 2 - parametertype%nqp = 2 + parametertype = simpleParameterType(1,2,2,0) + !parametertype%nodes_per_elem = 2 + !parametertype%nqp = 2 n = parametertype%nodes_per_elem call AllocAry(test_shape, parametertype%nodes_per_elem, parametertype%nqp, "test_shape", ErrStat, ErrMsg) @@ -96,10 +96,11 @@ subroutine test_BD_diffmtc_2node() deallocate(test_shape) deallocate(test_shapederivative) - deallocate(parametertype%QPtN) deallocate(gll_nodes) deallocate(baseline_shape) deallocate(baseline_shapederivative) + + call simpleparametertype_teardown(parametertype) end subroutine @@ -116,9 +117,9 @@ subroutine test_BD_diffmtc_5node() ! -------------------------------------------------------------------------- testname = "5-node element: evaluate shape/shapederivative at nodes and at three non-node locations" - parametertype = simpleParameterType() - parametertype%nodes_per_elem = 5 - parametertype%nqp = 8 ! testing the GLL nodes and three non-nodal locations (-0.8, 0.1, 0.4) + parametertype = simpleParameterType(1,5,8,0) + !parametertype%nodes_per_elem = 5 + !parametertype%nqp = 8 ! testing the GLL nodes and three non-nodal locations (-0.8, 0.1, 0.4) n = parametertype%nodes_per_elem call AllocAry(test_shape, parametertype%nodes_per_elem, parametertype%nqp, "test_shape", ErrStat, ErrMsg) @@ -154,10 +155,11 @@ subroutine test_BD_diffmtc_5node() deallocate(test_shape) deallocate(test_shapederivative) - deallocate(parametertype%QPtN) deallocate(gll_nodes) deallocate(baseline_shape) deallocate(baseline_shapederivative) + + call simpleparametertype_teardown(parametertype) end subroutine end module diff --git a/modules/beamdyn/tests/test_ExtractRelativeRotation.F90 b/modules/beamdyn/tests/test_ExtractRelativeRotation.F90 index 38bdea6eed..0d43c30c91 100644 --- a/modules/beamdyn/tests/test_ExtractRelativeRotation.F90 +++ b/modules/beamdyn/tests/test_ExtractRelativeRotation.F90 @@ -26,7 +26,7 @@ subroutine test_ExtractRelativeRotation() ! -------------------------------------------------------------------------- testname = "static simple beam under gravity:" - parametertype = simpleParameterType() + parametertype = simpleParameterType(1,16,16,0) call ExtractRelativeRotation(identity(), parametertype, rr, ErrStat, ErrMsg) diff --git a/modules/beamdyn/tests/test_tools.F90 b/modules/beamdyn/tests/test_tools.F90 index e0ea179533..45e4036ed0 100644 --- a/modules/beamdyn/tests/test_tools.F90 +++ b/modules/beamdyn/tests/test_tools.F90 @@ -91,17 +91,27 @@ function getGravityInZ() getGravityInZ = (/ 0.0, 0.0, -9.806 /) end function - type(BD_ParameterType) function simpleParameterType() RESULT(p) - + type(BD_ParameterType) function simpleParameterType(elem_total,nodes_per_elem,nqp,qp_indx_offset) RESULT(p) + + integer, intent(in ) :: elem_total + integer, intent(in ) :: nodes_per_elem + integer, intent(in ) :: nqp + integer, intent(in ) :: qp_indx_offset + integer :: i, j integer :: ErrStat character(1024) :: ErrMsg ! scalars - p%elem_total = 1 - p%nodes_per_elem = 16 - p%nqp = 16 - p%qp_indx_offset = 0 + !p%elem_total = 1 + !p%nodes_per_elem = 16 + !p%nqp = 16 + !p%qp_indx_offset = 0 + + p%elem_total = elem_total + p%nodes_per_elem = nodes_per_elem + p%nqp = nqp + p%qp_indx_offset = qp_indx_offset ! fixed size arrays p%Glb_crv = (/ 0.0, 0.0, 0.0 /) @@ -117,9 +127,11 @@ type(BD_ParameterType) function simpleParameterType() RESULT(p) call AllocAry(p%QPtw_Shp_Jac, p%nqp, p%nodes_per_elem, p%elem_total, 'QPtw_Shp_Jac', ErrStat, ErrMsg) call AllocAry(p%Shp, p%nodes_per_elem, p%nqp, 'Shp', ErrStat, ErrMsg) call AllocAry(p%ShpDer, p%nodes_per_elem, p%nqp, 'ShpDer', ErrStat, ErrMsg) + call AllocAry(p%QPtN, p%nqp, 'QPtN', ErrStat, ErrMsg) call AllocAry(p%QPtWeight, p%nqp, 'QPtWeightShp', ErrStat, ErrMsg) call AllocAry(p%QPtw_ShpDer, p%nqp, p%nodes_per_elem, 'QPtw_ShpDer', ErrStat, ErrMsg) - call AllocAry(p%Jacobian, p%nqp, p%nodes_per_elem, 'Jacobian', ErrStat, ErrMsg) + call AllocAry(p%Jacobian, p%nqp, p%elem_total, 'Jacobian', ErrStat, ErrMsg) + call AllocAry(p%uuN0, 3, p%nodes_per_elem, p%elem_total,'uuN0', ErrStat, ErrMsg) ! construct arrays p%qp%mmm = 1.0 @@ -132,7 +144,27 @@ type(BD_ParameterType) function simpleParameterType() RESULT(p) end do end function - + + subroutine SimpleParameterType_TearDown(p) + + TYPE(BD_ParameterType), intent(inout) :: p + + deallocate(p%qp%mmm) + deallocate(p%qp%mEta) + deallocate(p%Mass0_QP) + deallocate(p%QPtw_Shp_Shp_Jac) + deallocate(p%QPtw_ShpDer_ShpDer_Jac) + deallocate(p%QPtw_Shp_ShpDer) + deallocate(p%QPtw_Shp_Jac) + deallocate(p%Shp) + deallocate(p%ShpDer) + deallocate(p%QPtN) + deallocate(p%QPtWeight) + deallocate(p%QPtw_ShpDer) + deallocate(p%Jacobian) + deallocate(p%uuN0) + end subroutine + type(BD_MiscVarType) function simpleMiscVarType(nqp, nelem) RESULT(m) integer, intent(in) :: nqp, nelem From 7a995935929b25d8e4e2005a04d2a5b352ff7270 Mon Sep 17 00:00:00 2001 From: Mike Sprague Date: Fri, 9 Oct 2020 11:45:28 -0600 Subject: [PATCH 4/8] adding test: test_BD_MemberEta.f90 --- modules/beamdyn/tests/test_BD_MemberEta.F90 | 82 +++++++++++++++++++++ unit_tests/beamdyn/CMakeLists.txt | 1 + 2 files changed, 83 insertions(+) create mode 100644 modules/beamdyn/tests/test_BD_MemberEta.F90 diff --git a/modules/beamdyn/tests/test_BD_MemberEta.F90 b/modules/beamdyn/tests/test_BD_MemberEta.F90 new file mode 100644 index 0000000000..7dacd5944b --- /dev/null +++ b/modules/beamdyn/tests/test_BD_MemberEta.F90 @@ -0,0 +1,82 @@ +module test_BD_MemberEta + + ! tests routine that calculates the length of the beam's reference line + ! also finds element boundaries in eta coordinates + + use pFUnit_mod + use BeamDyn + use NWTC_Num + use test_tools + + implicit none + + integer(IntKi) :: nqp + integer(IntKi) :: member_total + real(BDKi) :: total_length + real(BDKi), allocatable :: baseline_jac(:,:) + real(BDKi), allocatable :: baseline_QPtW(:) + real(BDKi), allocatable :: baseline_member_eta(:) + real(BDKi), allocatable :: test_member_eta(:) + real(BDKi) :: baseline_total_length + real(BDKi) :: test_total_length + integer(IntKi) :: ErrStat + character :: ErrMsg + character(1024) :: testname + real(BDKi) :: tolerance + +contains + + @test + subroutine test_BD_MemberEta_5node() + + ! initialize NWTC_Num constants + call SetConstants() + + tolerance = 1e-14 + + ! -------------------------------------------------------------------------- + testname = "test_bd_member_eta_5node" + + ! this test problem is for a beam reference axis defined by 5 points + ! x = 0., 0.16237631096713473, 0.25, -0.30523345382427747, -1. + ! y = 0., 0.17578464768961147, 1., 2.4670724951675314, 3.5 + ! z = 0., 0.1481911137890286, 1.1875, 2.953849702537502, 4. + ! where the underlying forth-order polynomial is + ! fx[u_] = u - 2 u^3; + ! fy[u_] = u/2 + 3 u^2; + ! fz[u_] = 5 u^2 - u^4; + ! + ! exact (to mp) length is 5.627175237247959 + + nqp = 5 ! number of quadrature points + member_total = 1 ! 1 element + + call AllocAry(baseline_jac , nqp, member_total, 'Reference Jacobian', ErrStat, ErrMsg) + call AllocAry(baseline_QPtW, nqp, 'Reference QPtWeight', ErrStat, ErrMsg) + call AllocAry(baseline_member_eta, member_total, 'Reference member_eta', ErrStat, ErrMsg) + call AllocAry(test_member_eta, member_total, 'test member_eta', ErrStat, ErrMsg) + + ! 5-point Guass-Legendre quadrature; see https://pomax.github.io/bezierinfo/legendre-gauss.html + baseline_QPtW(1:nqp) = (/ 0.2369268850561891, 0.4786286704993665, 0.5688888888888889, 0.4786286704993665, 0.2369268850561891 /) + + ! assign baseline jacobian based; these values were calculated in separate mathematica script + baseline_jac(1:nqp,1) = (/ 0.6715870058501458, 1.509599209717604, 2.861380785564901, 4.097191592895223, 4.880926263217582 /) + + ! total length of beam calculated in mathematica + !baseline_total_length = 5.627175237247959 ! this is actual length based on mathematica + !baseline_total_length = 5.627202424388781 ! this is approximation with 7-point gauss quadrature; 0.0005% error + baseline_total_length = 5.626918236484061 ! this is approximation with 5-point gauss quadrature; 0.005% error + baseline_member_eta(1) = 1. ! just one element + + call BD_MemberEta(member_total, baseline_QPtW, baseline_jac, test_member_eta, test_total_length) + + @assertEqual(baseline_total_length, test_total_length, tolerance, testname) + @assertEqual(baseline_member_eta, test_member_eta, tolerance, testname) + + deallocate(baseline_Jac) + deallocate(baseline_QPtW) + deallocate(baseline_member_eta) + deallocate(test_member_eta) + + end subroutine +end module diff --git a/unit_tests/beamdyn/CMakeLists.txt b/unit_tests/beamdyn/CMakeLists.txt index f998d29c55..aac7fdb66a 100644 --- a/unit_tests/beamdyn/CMakeLists.txt +++ b/unit_tests/beamdyn/CMakeLists.txt @@ -45,6 +45,7 @@ set(testlist test_BD_GaussPointWeight test_BD_diffmtc test_BD_InitShpDerJaco + test_BD_MemberEta ) foreach(test ${testlist}) set(test_dependency pfunit ${source_modulesdirectory}/${module_directory}/tests/${test}.F90) From 20c93abedd009e2c0af4f15a7ccc2e7c01247ea6 Mon Sep 17 00:00:00 2001 From: Mike Sprague Date: Fri, 9 Oct 2020 12:12:40 -0600 Subject: [PATCH 5/8] adding more comments to test_BD_MemberEta.F90 --- modules/beamdyn/tests/test_BD_MemberEta.F90 | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/modules/beamdyn/tests/test_BD_MemberEta.F90 b/modules/beamdyn/tests/test_BD_MemberEta.F90 index 7dacd5944b..08df56dbfd 100644 --- a/modules/beamdyn/tests/test_BD_MemberEta.F90 +++ b/modules/beamdyn/tests/test_BD_MemberEta.F90 @@ -47,6 +47,8 @@ subroutine test_BD_MemberEta_5node() ! fz[u_] = 5 u^2 - u^4; ! ! exact (to mp) length is 5.627175237247959 + ! The Jacobian values below are calculate based on a 5-node Legendre Spectral Element + ! While we give baseline Jacobian here, we check "test_BD_InitShpDerJaco.F90" elsewhere nqp = 5 ! number of quadrature points member_total = 1 ! 1 element @@ -62,11 +64,14 @@ subroutine test_BD_MemberEta_5node() ! assign baseline jacobian based; these values were calculated in separate mathematica script baseline_jac(1:nqp,1) = (/ 0.6715870058501458, 1.509599209717604, 2.861380785564901, 4.097191592895223, 4.880926263217582 /) - ! total length of beam calculated in mathematica + ! total length of beam calculated in mathematica; note that for this curved beam, GL quadrature is APPROXIMATE, but + ! converged rapidly; TR quadrature is not nearly as good; commented-out values are useful for understanding quadrature performance + ! for curved beams. !baseline_total_length = 5.627175237247959 ! this is actual length based on mathematica + !baseline_total_length = 5.634413547964786 ! this is approximation with 9-point Trapezoidal quadrature; 0.13% error !baseline_total_length = 5.627202424388781 ! this is approximation with 7-point gauss quadrature; 0.0005% error - baseline_total_length = 5.626918236484061 ! this is approximation with 5-point gauss quadrature; 0.005% error - baseline_member_eta(1) = 1. ! just one element + baseline_total_length = 5.626918236484061 ! this is approximation with 5-point gauss quadrature; 0.005% error (tested here) + baseline_member_eta(1) = 1. ! just one element; so member_length / total_length = 1 call BD_MemberEta(member_total, baseline_QPtW, baseline_jac, test_member_eta, test_total_length) From 77e0c5c80431a17962bafc51aa62498068471a88 Mon Sep 17 00:00:00 2001 From: Mike Sprague Date: Fri, 9 Oct 2020 16:31:36 -0600 Subject: [PATCH 6/8] adding test for BD_TrapezoidalPointWeight.F90; modifed target function in BeamDyn to simplify call --- modules/beamdyn/src/BeamDyn.f90 | 2 +- modules/beamdyn/src/BeamDyn_Subs.f90 | 19 ++- .../beamdyn/tests/test_BD_DistrLoadCopy.F90 | 2 +- .../beamdyn/tests/test_BD_GravityForce.F90 | 2 +- .../beamdyn/tests/test_BD_InitShpDerJaco.F90 | 2 +- .../beamdyn/tests/test_BD_QPData_mEta_rho.F90 | 2 +- .../tests/test_BD_TrapezoidalPointWeight.F90 | 113 ++++++++++++++++++ modules/beamdyn/tests/test_BD_diffmtc.F90 | 4 +- .../tests/test_ExtractRelativeRotation.F90 | 2 +- modules/beamdyn/tests/test_tools.F90 | 6 +- unit_tests/beamdyn/CMakeLists.txt | 1 + 11 files changed, 135 insertions(+), 20 deletions(-) create mode 100644 modules/beamdyn/tests/test_BD_TrapezoidalPointWeight.F90 diff --git a/modules/beamdyn/src/BeamDyn.f90 b/modules/beamdyn/src/BeamDyn.f90 index 7af5c177a2..975f69554c 100644 --- a/modules/beamdyn/src/BeamDyn.f90 +++ b/modules/beamdyn/src/BeamDyn.f90 @@ -152,7 +152,7 @@ SUBROUTINE BD_Init( InitInp, u, p, x, xd, z, OtherState, y, MiscVar, Interval, I ELSEIF(p%quadrature .EQ. TRAP_QUADRATURE) THEN - CALL BD_TrapezoidalPointWeight(p, InputFileData) ! computes p%QPtN and p%QPtWeight + CALL BD_TrapezoidalPointWeight(p, InputFileData%InpBl%station_eta, InputFileData%InpBl%station_total) ! computes p%QPtN and p%QPtWeight ENDIF diff --git a/modules/beamdyn/src/BeamDyn_Subs.f90 b/modules/beamdyn/src/BeamDyn_Subs.f90 index 9641704a9e..a76f9661a4 100644 --- a/modules/beamdyn/src/BeamDyn_Subs.f90 +++ b/modules/beamdyn/src/BeamDyn_Subs.f90 @@ -580,10 +580,11 @@ SUBROUTINE BD_GaussPointWeight(n, x, w, ErrStat, ErrMsg) END SUBROUTINE BD_GaussPointWeight !----------------------------------------------------------------------------------------------------------------------------------- ! This subroutine computes trapezoidal quadrature points and weights, p%QPtN and p%QPtWeight -SUBROUTINE BD_TrapezoidalPointWeight(p, InputFileData) +SUBROUTINE BD_TrapezoidalPointWeight(p, station_eta, station_total) TYPE(BD_ParameterType),INTENT(INOUT):: p !< BeamDyn parameters - TYPE(BD_InputFile), INTENT(IN ):: InputFileData !< BeamDyn input-file data + Integer(IntKi),INTENT(IN ) :: station_total + REAL(BDKi),INTENT(IN ) :: station_eta(:) ! local variables REAL(BDKi) :: denom ! denominator for quadrature weight computations @@ -593,20 +594,20 @@ SUBROUTINE BD_TrapezoidalPointWeight(p, InputFileData) INTEGER(IntKi) :: id1, j !bjj: this assumes there is only one member - + + ! compute the trapezoidal quadrature points, p%QPtN, and scale to range [-1,1]: ! If there is refinement, this will add new points between the specified ones. If p%refine == 1, can skip this. - p%QPtN(1) = InputFileData%InpBl%station_eta(1) + p%QPtN(1) = station_eta(1) DO j = 2,p%nqp indx = 1+(j-2_IntKi)/p%refine ! note use of integer math here --> (J-2)/p%refine may not be integer. - p%QPtN(j) = InputFileData%InpBl%station_eta(indx) + & - ((InputFileData%InpBl%station_eta(indx+1) - InputFileData%InpBl%station_eta(indx))/p%refine) * (MOD(j-2,p%refine) + 1) + p%QPtN(j) = station_eta(indx) + & + ((station_eta(indx+1) - station_eta(indx))/p%refine) * (MOD(j-2,p%refine) + 1) ENDDO p%QPtN = 2.0_BDKi*p%QPtN - 1.0_BDKi ! rescale range from [0, 1] to [-1,1] - ! compute the trapezoidal quadrature weights, p%QPtWeight: - id1 = InputFileData%InpBl%station_total + id1 = station_total temp_id0 = (id0 - 1)*p%refine + 1 ! Starting index in QPtN --> always going to be 1 temp_id1 = (id1 - 1)*p%refine + 1 ! ending index in QPtN --> will be size(p%QPtN) denom = p%QPtN(temp_id1) - p%QPtN(temp_id0) ! This is the range of QPtN --> for single member, is always == 2 @@ -617,8 +618,6 @@ SUBROUTINE BD_TrapezoidalPointWeight(p, InputFileData) ENDDO p%QPtWeight(p%nqp) = (p%QPtN(temp_id1 ) - p%QPtN(temp_id1-1 ))/denom - - END SUBROUTINE BD_TrapezoidalPointWeight !----------------------------------------------------------------------------------------------------------------------------------- diff --git a/modules/beamdyn/tests/test_BD_DistrLoadCopy.F90 b/modules/beamdyn/tests/test_BD_DistrLoadCopy.F90 index 88e5ae87f6..dbc68fbd4f 100644 --- a/modules/beamdyn/tests/test_BD_DistrLoadCopy.F90 +++ b/modules/beamdyn/tests/test_BD_DistrLoadCopy.F90 @@ -29,7 +29,7 @@ subroutine test_BD_DistrLoadCopy() testname = "static simple beam under gravity:" ! build the parametertype, inputtype, and miscvartype - parametertype = simpleParameterType(1,16,16,0) + parametertype = simpleParameterType(1,16,16,0,1) miscvartype = simpleMiscVarType(parametertype%nqp, parametertype%elem_total) inputtype = simpleInputType(parametertype%nqp, parametertype%elem_total) diff --git a/modules/beamdyn/tests/test_BD_GravityForce.F90 b/modules/beamdyn/tests/test_BD_GravityForce.F90 index aa99473122..d5b54a762b 100644 --- a/modules/beamdyn/tests/test_BD_GravityForce.F90 +++ b/modules/beamdyn/tests/test_BD_GravityForce.F90 @@ -30,7 +30,7 @@ subroutine test_BD_GravityForce() baseline(4:6) = (/ 0.0, 0.0, 0.0 /) ! allocate and build the custom types - parametertype = simpleParameterType(1,16,16,0) + parametertype = simpleParameterType(1,16,16,0,1) miscvartype = simpleMiscVarType(parametertype%nqp, parametertype%elem_total) gravity = getGravityInZ() diff --git a/modules/beamdyn/tests/test_BD_InitShpDerJaco.F90 b/modules/beamdyn/tests/test_BD_InitShpDerJaco.F90 index a1e88cb038..f2f8877853 100644 --- a/modules/beamdyn/tests/test_BD_InitShpDerJaco.F90 +++ b/modules/beamdyn/tests/test_BD_InitShpDerJaco.F90 @@ -34,7 +34,7 @@ subroutine test_BD_InitShpDerJaco_5node() ! -------------------------------------------------------------------------- testname = "5 node, 1 element, curved:" - p = simpleparametertype(1,5,5,0) + p = simpleparametertype(1,5,5,0,1) !p%elem_total = 1 !p%nodes_per_elem = 5 !p%nqp = 5 ! Let's use Gauss_Legendre Quadrature, which should be exact for intended polynomial test case diff --git a/modules/beamdyn/tests/test_BD_QPData_mEta_rho.F90 b/modules/beamdyn/tests/test_BD_QPData_mEta_rho.F90 index 3382904b2d..58012d1ee5 100644 --- a/modules/beamdyn/tests/test_BD_QPData_mEta_rho.F90 +++ b/modules/beamdyn/tests/test_BD_QPData_mEta_rho.F90 @@ -34,7 +34,7 @@ subroutine test_BD_QPData_mEta_rho() baselineRR0mEta = (/ 0.0, 0.0, 0.0 /) ! allocate and build the custom input types - parametertype = simpleParameterType(1,16,16,0) + parametertype = simpleParameterType(1,16,16,0,1) miscvartype = simpleMiscVarType(parametertype%nqp, parametertype%elem_total) ! allocate the results diff --git a/modules/beamdyn/tests/test_BD_TrapezoidalPointWeight.F90 b/modules/beamdyn/tests/test_BD_TrapezoidalPointWeight.F90 new file mode 100644 index 0000000000..64791336f6 --- /dev/null +++ b/modules/beamdyn/tests/test_BD_TrapezoidalPointWeight.F90 @@ -0,0 +1,113 @@ +module test_BD_TrapezoidalPointWeight + + use pFUnit_mod + use BeamDyn + use NWTC_Num + use test_tools + + implicit none + + type(BD_ParameterType) :: p + + integer(IntKi) :: nqp + integer(IntKi) :: refine + integer(IntKi) :: station_total + real(BDKi), allocatable :: station_eta(:) + real(BDKi), allocatable :: baseline_QPtN(:) + real(BDKi), allocatable :: baseline_QPtW(:) + + integer(IntKi) :: ErrStat + character :: ErrMsg + character(1024) :: testname + real(BDKi) :: tolerance + +contains + + @test + subroutine test_BD_TrapezoidalPointWeight_2station() + + ! initialize NWTC_Num constants + call SetConstants() + + tolerance = 1e-14 + + ! -------------------------------------------------------------------------- + testname = "test_BD_TrapezoidalPointWeight_2station_8refine" + + station_total = 2 + + call AllocAry(station_eta, station_total, "station_eta", ErrStat, ErrMsg) + + ! simple case where we have enpoints only; typical with constant cross sections + station_eta(1:station_total) = (/0., 1./) + + refine = 8 + + nqp = (station_total - 1)*refine + 1 + + p=simpleParameterType(1,1,nqp,0,refine) + + call AllocAry(baseline_QPtN, nqp, "baseline_QPtN", ErrStat, ErrMsg) + call AllocAry(baseline_QPtW, nqp, "baseline_QPtW", ErrStat, ErrMsg) + + baseline_QPtN(1:nqp) = (/ -1., -0.75, -0.5, -0.25, 0., 0.25, 0.5, 0.75, 1. /) + baseline_QPtW(1:nqp) = (/ 0.125, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.125 /) + + call BD_TrapezoidalPointWeight(p, station_eta, station_total) + + @assertEqual(baseline_QPtN, p%QPtN, tolerance, testname) + @assertEqual(baseline_QPtW, p%QPtWeight, tolerance, testname) + + deallocate(station_eta) + deallocate(baseline_QPtN) + deallocate(baseline_QPtW) + + call simpleparametertype_teardown(p) + + end subroutine + + @test + subroutine test_BD_TrapezoidalPointWeight_3station() + + ! initialize NWTC_Num constants + call SetConstants() + + tolerance = 1e-14 + + ! -------------------------------------------------------------------------- + testname = "test_BD_TrapezoidalPointWeight_3station_2refine" + + ! provide three stations, unequally distributed + ! refine by factor of two + + station_total = 3 + + call AllocAry(station_eta, station_total, "station_eta", ErrStat, ErrMsg) + + station_eta(1:station_total) = (/0., 0.25, 1./) + + refine = 2 + + nqp = (station_total - 1)*refine + 1 + + p=simpleParameterType(1,1,nqp,0,refine) + + call AllocAry(baseline_QPtN, nqp, "baseline_QPtN", ErrStat, ErrMsg) + call AllocAry(baseline_QPtW, nqp, "baseline_QPtW", ErrStat, ErrMsg) + + baseline_QPtN(1:nqp) = (/ -1., -0.75, -0.5, 0.25, 1. /) + baseline_QPtW(1:nqp) = (/ 0.125, 0.25, 0.125+0.5*0.75, 0.75, 0.5*0.75/) + + call BD_TrapezoidalPointWeight(p, station_eta, station_total) + + @assertEqual(baseline_QPtN, p%QPtN, tolerance, testname) + @assertEqual(baseline_QPtW, p%QPtWeight, tolerance, testname) + + deallocate(station_eta) + deallocate(baseline_QPtN) + deallocate(baseline_QPtW) + + call simpleparametertype_teardown(p) + + end subroutine +end module diff --git a/modules/beamdyn/tests/test_BD_diffmtc.F90 b/modules/beamdyn/tests/test_BD_diffmtc.F90 index 84093e5423..9a1dfa88f2 100644 --- a/modules/beamdyn/tests/test_BD_diffmtc.F90 +++ b/modules/beamdyn/tests/test_BD_diffmtc.F90 @@ -68,7 +68,7 @@ subroutine test_BD_diffmtc_2node() ! shpder(1,:) = -0.5, -0.5 ! shpder(2,:) = 0.5, 0.5 - parametertype = simpleParameterType(1,2,2,0) + parametertype = simpleParameterType(1,2,2,0,1) !parametertype%nodes_per_elem = 2 !parametertype%nqp = 2 n = parametertype%nodes_per_elem @@ -117,7 +117,7 @@ subroutine test_BD_diffmtc_5node() ! -------------------------------------------------------------------------- testname = "5-node element: evaluate shape/shapederivative at nodes and at three non-node locations" - parametertype = simpleParameterType(1,5,8,0) + parametertype = simpleParameterType(1,5,8,0,1) !parametertype%nodes_per_elem = 5 !parametertype%nqp = 8 ! testing the GLL nodes and three non-nodal locations (-0.8, 0.1, 0.4) n = parametertype%nodes_per_elem diff --git a/modules/beamdyn/tests/test_ExtractRelativeRotation.F90 b/modules/beamdyn/tests/test_ExtractRelativeRotation.F90 index 0d43c30c91..25c1c01dcf 100644 --- a/modules/beamdyn/tests/test_ExtractRelativeRotation.F90 +++ b/modules/beamdyn/tests/test_ExtractRelativeRotation.F90 @@ -26,7 +26,7 @@ subroutine test_ExtractRelativeRotation() ! -------------------------------------------------------------------------- testname = "static simple beam under gravity:" - parametertype = simpleParameterType(1,16,16,0) + parametertype = simpleParameterType(1,16,16,0,0) call ExtractRelativeRotation(identity(), parametertype, rr, ErrStat, ErrMsg) diff --git a/modules/beamdyn/tests/test_tools.F90 b/modules/beamdyn/tests/test_tools.F90 index 45e4036ed0..911709b35f 100644 --- a/modules/beamdyn/tests/test_tools.F90 +++ b/modules/beamdyn/tests/test_tools.F90 @@ -91,12 +91,13 @@ function getGravityInZ() getGravityInZ = (/ 0.0, 0.0, -9.806 /) end function - type(BD_ParameterType) function simpleParameterType(elem_total,nodes_per_elem,nqp,qp_indx_offset) RESULT(p) + type(BD_ParameterType) function simpleParameterType(elem_total,nodes_per_elem,nqp,qp_indx_offset,refine) RESULT(p) integer, intent(in ) :: elem_total integer, intent(in ) :: nodes_per_elem integer, intent(in ) :: nqp integer, intent(in ) :: qp_indx_offset + integer, intent(in ) :: refine integer :: i, j integer :: ErrStat @@ -112,6 +113,7 @@ type(BD_ParameterType) function simpleParameterType(elem_total,nodes_per_elem,nq p%nodes_per_elem = nodes_per_elem p%nqp = nqp p%qp_indx_offset = qp_indx_offset + p%refine = refine ! fixed size arrays p%Glb_crv = (/ 0.0, 0.0, 0.0 /) @@ -128,7 +130,7 @@ type(BD_ParameterType) function simpleParameterType(elem_total,nodes_per_elem,nq call AllocAry(p%Shp, p%nodes_per_elem, p%nqp, 'Shp', ErrStat, ErrMsg) call AllocAry(p%ShpDer, p%nodes_per_elem, p%nqp, 'ShpDer', ErrStat, ErrMsg) call AllocAry(p%QPtN, p%nqp, 'QPtN', ErrStat, ErrMsg) - call AllocAry(p%QPtWeight, p%nqp, 'QPtWeightShp', ErrStat, ErrMsg) + call AllocAry(p%QPtWeight, p%nqp, 'QPtWeight', ErrStat, ErrMsg) call AllocAry(p%QPtw_ShpDer, p%nqp, p%nodes_per_elem, 'QPtw_ShpDer', ErrStat, ErrMsg) call AllocAry(p%Jacobian, p%nqp, p%elem_total, 'Jacobian', ErrStat, ErrMsg) call AllocAry(p%uuN0, 3, p%nodes_per_elem, p%elem_total,'uuN0', ErrStat, ErrMsg) diff --git a/unit_tests/beamdyn/CMakeLists.txt b/unit_tests/beamdyn/CMakeLists.txt index aac7fdb66a..6adf133c28 100644 --- a/unit_tests/beamdyn/CMakeLists.txt +++ b/unit_tests/beamdyn/CMakeLists.txt @@ -46,6 +46,7 @@ set(testlist test_BD_diffmtc test_BD_InitShpDerJaco test_BD_MemberEta + test_BD_TrapezoidalPointWeight ) foreach(test ${testlist}) set(test_dependency pfunit ${source_modulesdirectory}/${module_directory}/tests/${test}.F90) From 4c5ce546427af725e8181943811bdc49d2899f6d Mon Sep 17 00:00:00 2001 From: Rafael M Mudafort Date: Wed, 14 Oct 2020 10:37:48 -0500 Subject: [PATCH 7/8] Remove commented and obsolete code --- .../beamdyn/tests/test_BD_InitShpDerJaco.F90 | 4 +--- modules/beamdyn/tests/test_BD_diffmtc.F90 | 2 -- modules/beamdyn/tests/test_tools.F90 | 18 ++++++------------ 3 files changed, 7 insertions(+), 17 deletions(-) diff --git a/modules/beamdyn/tests/test_BD_InitShpDerJaco.F90 b/modules/beamdyn/tests/test_BD_InitShpDerJaco.F90 index f2f8877853..78c0080030 100644 --- a/modules/beamdyn/tests/test_BD_InitShpDerJaco.F90 +++ b/modules/beamdyn/tests/test_BD_InitShpDerJaco.F90 @@ -34,10 +34,8 @@ subroutine test_BD_InitShpDerJaco_5node() ! -------------------------------------------------------------------------- testname = "5 node, 1 element, curved:" + ! Let's use Gauss_Legendre Quadrature, which should be exact for intended polynomial test case p = simpleparametertype(1,5,5,0,1) - !p%elem_total = 1 - !p%nodes_per_elem = 5 - !p%nqp = 5 ! Let's use Gauss_Legendre Quadrature, which should be exact for intended polynomial test case ! Allocate memory for baseline results call AllocAry(baseline_Shp , p%nodes_per_elem, p%nqp, 'Reference Shp' , ErrStat, ErrMsg) diff --git a/modules/beamdyn/tests/test_BD_diffmtc.F90 b/modules/beamdyn/tests/test_BD_diffmtc.F90 index 9a1dfa88f2..ccf6584f70 100644 --- a/modules/beamdyn/tests/test_BD_diffmtc.F90 +++ b/modules/beamdyn/tests/test_BD_diffmtc.F90 @@ -69,8 +69,6 @@ subroutine test_BD_diffmtc_2node() ! shpder(2,:) = 0.5, 0.5 parametertype = simpleParameterType(1,2,2,0,1) - !parametertype%nodes_per_elem = 2 - !parametertype%nqp = 2 n = parametertype%nodes_per_elem call AllocAry(test_shape, parametertype%nodes_per_elem, parametertype%nqp, "test_shape", ErrStat, ErrMsg) diff --git a/modules/beamdyn/tests/test_tools.F90 b/modules/beamdyn/tests/test_tools.F90 index 911709b35f..1a5f28f9a5 100644 --- a/modules/beamdyn/tests/test_tools.F90 +++ b/modules/beamdyn/tests/test_tools.F90 @@ -91,34 +91,28 @@ function getGravityInZ() getGravityInZ = (/ 0.0, 0.0, -9.806 /) end function - type(BD_ParameterType) function simpleParameterType(elem_total,nodes_per_elem,nqp,qp_indx_offset,refine) RESULT(p) + type(BD_ParameterType) function simpleParameterType(elem_total, nodes_per_elem, nqp, qp_indx_offset, refine) RESULT(p) integer, intent(in ) :: elem_total integer, intent(in ) :: nodes_per_elem integer, intent(in ) :: nqp integer, intent(in ) :: qp_indx_offset integer, intent(in ) :: refine - + integer :: i, j integer :: ErrStat character(1024) :: ErrMsg - - ! scalars - !p%elem_total = 1 - !p%nodes_per_elem = 16 - !p%nqp = 16 - !p%qp_indx_offset = 0 p%elem_total = elem_total p%nodes_per_elem = nodes_per_elem p%nqp = nqp p%qp_indx_offset = qp_indx_offset p%refine = refine - + ! fixed size arrays p%Glb_crv = (/ 0.0, 0.0, 0.0 /) p%GlbRot = identity() - + ! allocate arrays call AllocAry(p%qp%mmm, p%nqp, p%elem_total, 'qp_mmm', ErrStat, ErrMsg) call AllocAry(p%qp%mEta, 3, p%nqp, p%elem_total, 'qp_RR0mEta', ErrStat, ErrMsg) @@ -137,14 +131,14 @@ type(BD_ParameterType) function simpleParameterType(elem_total,nodes_per_elem,nq ! construct arrays p%qp%mmm = 1.0 - + do j=1, p%elem_total do i=1, p%nqp p%qp%mEta(:,i,j) = (/ 0.0, 0.0, 0.0 /) p%Mass0_QP(:,:,(i-1)*p%elem_total+j) = getMassMatrix() end do end do - + end function subroutine SimpleParameterType_TearDown(p) From 2974f82bc67aeee5fb12462a004452de21a5a08b Mon Sep 17 00:00:00 2001 From: Rafael M Mudafort Date: Wed, 14 Oct 2020 11:04:12 -0500 Subject: [PATCH 8/8] Use the BD Types deallocate method --- .../beamdyn/tests/test_BD_DistrLoadCopy.F90 | 2 +- .../beamdyn/tests/test_BD_GravityForce.F90 | 2 +- .../beamdyn/tests/test_BD_InitShpDerJaco.F90 | 2 +- .../beamdyn/tests/test_BD_QPData_mEta_rho.F90 | 2 +- .../tests/test_BD_TrapezoidalPointWeight.F90 | 4 ++-- modules/beamdyn/tests/test_BD_diffmtc.F90 | 4 ++-- modules/beamdyn/tests/test_tools.F90 | 20 ------------------- 7 files changed, 8 insertions(+), 28 deletions(-) diff --git a/modules/beamdyn/tests/test_BD_DistrLoadCopy.F90 b/modules/beamdyn/tests/test_BD_DistrLoadCopy.F90 index dbc68fbd4f..bfa6037513 100644 --- a/modules/beamdyn/tests/test_BD_DistrLoadCopy.F90 +++ b/modules/beamdyn/tests/test_BD_DistrLoadCopy.F90 @@ -42,5 +42,5 @@ subroutine test_BD_DistrLoadCopy() end do end do - call simpleparametertype_teardown(parametertype) + call BD_DestroyParam(parametertype, ErrStat, ErrMsg) end subroutine diff --git a/modules/beamdyn/tests/test_BD_GravityForce.F90 b/modules/beamdyn/tests/test_BD_GravityForce.F90 index d5b54a762b..6b7b4bc4dd 100644 --- a/modules/beamdyn/tests/test_BD_GravityForce.F90 +++ b/modules/beamdyn/tests/test_BD_GravityForce.F90 @@ -41,6 +41,6 @@ subroutine test_BD_GravityForce() ! test the values @assertEqual(baseline, miscvartype%qp%Fg(:,1,1), tolerance, testname) - call simpleparametertype_teardown(parametertype) + call BD_DestroyParam(parametertype, ErrStat, ErrMsg) end subroutine diff --git a/modules/beamdyn/tests/test_BD_InitShpDerJaco.F90 b/modules/beamdyn/tests/test_BD_InitShpDerJaco.F90 index 78c0080030..13e39f733b 100644 --- a/modules/beamdyn/tests/test_BD_InitShpDerJaco.F90 +++ b/modules/beamdyn/tests/test_BD_InitShpDerJaco.F90 @@ -167,7 +167,7 @@ subroutine test_BD_InitShpDerJaco_5node() deallocate(baseline_QPtw_Shp_Jac) deallocate(baseline_QPtw_ShpDer) - call SimpleParameterType_TearDown(p) + call BD_DestroyParam(p, ErrStat, ErrMsg) end subroutine end module diff --git a/modules/beamdyn/tests/test_BD_QPData_mEta_rho.F90 b/modules/beamdyn/tests/test_BD_QPData_mEta_rho.F90 index 58012d1ee5..5cea831001 100644 --- a/modules/beamdyn/tests/test_BD_QPData_mEta_rho.F90 +++ b/modules/beamdyn/tests/test_BD_QPData_mEta_rho.F90 @@ -46,5 +46,5 @@ subroutine test_BD_QPData_mEta_rho() @assertEqual(baselineRR0mEta, miscvartype%qp%RR0mEta(:,i,j), tolerance, testname) end do end do - call simpleParameterType_teardown(parametertype) + call BD_DestroyParam(parametertype, ErrStat, ErrMsg) end subroutine diff --git a/modules/beamdyn/tests/test_BD_TrapezoidalPointWeight.F90 b/modules/beamdyn/tests/test_BD_TrapezoidalPointWeight.F90 index 64791336f6..0ca075ad1a 100644 --- a/modules/beamdyn/tests/test_BD_TrapezoidalPointWeight.F90 +++ b/modules/beamdyn/tests/test_BD_TrapezoidalPointWeight.F90 @@ -62,7 +62,7 @@ subroutine test_BD_TrapezoidalPointWeight_2station() deallocate(baseline_QPtN) deallocate(baseline_QPtW) - call simpleparametertype_teardown(p) + call BD_DestroyParam(p, ErrStat, ErrMsg) end subroutine @@ -107,7 +107,7 @@ subroutine test_BD_TrapezoidalPointWeight_3station() deallocate(baseline_QPtN) deallocate(baseline_QPtW) - call simpleparametertype_teardown(p) + call BD_DestroyParam(p, ErrStat, ErrMsg) end subroutine end module diff --git a/modules/beamdyn/tests/test_BD_diffmtc.F90 b/modules/beamdyn/tests/test_BD_diffmtc.F90 index ccf6584f70..6d3319d92e 100644 --- a/modules/beamdyn/tests/test_BD_diffmtc.F90 +++ b/modules/beamdyn/tests/test_BD_diffmtc.F90 @@ -98,7 +98,7 @@ subroutine test_BD_diffmtc_2node() deallocate(baseline_shape) deallocate(baseline_shapederivative) - call simpleparametertype_teardown(parametertype) + call BD_DestroyParam(parametertype, ErrStat, ErrMsg) end subroutine @@ -157,7 +157,7 @@ subroutine test_BD_diffmtc_5node() deallocate(baseline_shape) deallocate(baseline_shapederivative) - call simpleparametertype_teardown(parametertype) + call BD_DestroyParam(parametertype, ErrStat, ErrMsg) end subroutine end module diff --git a/modules/beamdyn/tests/test_tools.F90 b/modules/beamdyn/tests/test_tools.F90 index 1a5f28f9a5..0b63f9fd4b 100644 --- a/modules/beamdyn/tests/test_tools.F90 +++ b/modules/beamdyn/tests/test_tools.F90 @@ -141,26 +141,6 @@ type(BD_ParameterType) function simpleParameterType(elem_total, nodes_per_elem, end function - subroutine SimpleParameterType_TearDown(p) - - TYPE(BD_ParameterType), intent(inout) :: p - - deallocate(p%qp%mmm) - deallocate(p%qp%mEta) - deallocate(p%Mass0_QP) - deallocate(p%QPtw_Shp_Shp_Jac) - deallocate(p%QPtw_ShpDer_ShpDer_Jac) - deallocate(p%QPtw_Shp_ShpDer) - deallocate(p%QPtw_Shp_Jac) - deallocate(p%Shp) - deallocate(p%ShpDer) - deallocate(p%QPtN) - deallocate(p%QPtWeight) - deallocate(p%QPtw_ShpDer) - deallocate(p%Jacobian) - deallocate(p%uuN0) - end subroutine - type(BD_MiscVarType) function simpleMiscVarType(nqp, nelem) RESULT(m) integer, intent(in) :: nqp, nelem