From 5a73dfca557a3385f8c740eb3252e5ef5bf8b838 Mon Sep 17 00:00:00 2001 From: Surya K S Date: Fri, 9 Aug 2024 15:58:27 +0200 Subject: [PATCH 1/3] Add all Dunavant quadrature rules for the triangle Currently, the first 7 rules of the triangle quadrature are Dunavant's rules. The higher order rules stem from somewhere else (TODO: Add reference when figured out) In addition to the existing mixed rule, we added Dunavant until order 20. The default is still the mixed rule. By wrapping the order in the type TriangleQuadDunavant, the user can select the Dunavant rule. TODO: There are more and better rules out there for integrating on a triangle (some Dunavant rules have negative weights) See, for example, the QuadPy-Package. Probably, we want a similar package for Julia in the future. --- src/CompScienceMeshes.jl | 2 +- src/quadrature/triangle/TriangleDunavant.jl | 2137 +++++++++++++++++ .../{ => triangle}/TriangleGauss.jl | 212 +- .../triangle/TriangleQuadratures.jl | 43 + test/runtests.jl | 1 + test/test_trdunavant.jl | 91 + test/test_trgauss.jl | 2 +- 7 files changed, 2375 insertions(+), 113 deletions(-) create mode 100644 src/quadrature/triangle/TriangleDunavant.jl rename src/quadrature/{ => triangle}/TriangleGauss.jl (97%) create mode 100644 src/quadrature/triangle/TriangleQuadratures.jl create mode 100644 test/test_trdunavant.jl diff --git a/src/CompScienceMeshes.jl b/src/CompScienceMeshes.jl index 56d66ef..e2d9715 100644 --- a/src/CompScienceMeshes.jl +++ b/src/CompScienceMeshes.jl @@ -95,7 +95,7 @@ include("utils/circumctr.jl") # quadrature rules for segements, triangles, and squares include("quadrature/SegmentGauss.jl") -include("quadrature/TriangleGauss.jl") +include("quadrature/triangle/TriangleQuadratures.jl") include("quadrature/SquareGauss.jl") # mesh component diff --git a/src/quadrature/triangle/TriangleDunavant.jl b/src/quadrature/triangle/TriangleDunavant.jl new file mode 100644 index 0000000..71244d6 --- /dev/null +++ b/src/quadrature/triangle/TriangleDunavant.jl @@ -0,0 +1,2137 @@ +# Cubature rules based on the paper +# Dunavant, D. A. “High Degree Efficient Symmetrical Gaussian Quadrature +# Rules for the Triangle.” International Journal for Numerical Methods +# in Engineering 21, no. 6 (1985): 1129–48. https://doi.org/10.1002/nme.1620210612. + +const trianglequadDunavant1a = @SVector [ + 0.333333333333333, +] + +const trianglequadDunavant1b = @SVector [ + 0.333333333333333, +] + +const trianglequadDunavant1w = @SVector [ + 1.000000000000000, +] + +const trianglequadDunavant3a = @SVector [ + 0.666666666666667, + 0.166666666666667, + 0.166666666666667, +] + +const trianglequadDunavant3b = @SVector [ + 0.166666666666667, + 0.166666666666667, + 0.666666666666667, +] + +const trianglequadDunavant3w = @SVector [ + 0.333333333333333, + 0.333333333333333, + 0.333333333333333, +] + +const trianglequadDunavant4a = @SVector [ + 0.333333333333333, + 0.600000000000000, + 0.200000000000000, + 0.200000000000000, +] + +const trianglequadDunavant4b = @SVector [ + 0.333333333333333, + 0.200000000000000, + 0.200000000000000, + 0.600000000000000, +] + +const trianglequadDunavant4w = @SVector [ + -0.562500000000000, + 0.520833333333333, + 0.520833333333333, + 0.520833333333333, +] + +const trianglequadDunavant6a = @SVector [ + 0.108103018168070, + 0.445948490915965, + 0.445948490915965, + 0.816847572980459, + 0.091576213509771, + 0.091576213509771, +] + +const trianglequadDunavant6b = @SVector [ + 0.445948490915965, + 0.445948490915965, + 0.108103018168070, + 0.091576213509771, + 0.091576213509771, + 0.816847572980459, +] + +const trianglequadDunavant6w = @SVector [ + 0.223381589678011, + 0.223381589678011, + 0.223381589678011, + 0.109951743655322, + 0.109951743655322, + 0.109951743655322, +] + +const trianglequadDunavant7a = @SVector [ + 0.333333333333333, + 0.059715871789770, + 0.470142064105115, + 0.470142064105115, + 0.797426985353087, + 0.101286507323456, + 0.101286507323456, +] + +const trianglequadDunavant7b = @SVector [ + 0.333333333333333, + 0.470142064105115, + 0.470142064105115, + 0.059715871789770, + 0.101286507323456, + 0.101286507323456, + 0.797426985353087, +] + +const trianglequadDunavant7w = @SVector [ + 0.225000000000000, + 0.132394152788506, + 0.132394152788506, + 0.132394152788506, + 0.125939180544827, + 0.125939180544827, + 0.125939180544827, +] + +const trianglequadDunavant12a = @SVector [ + 0.501426509658179, + 0.249286745170910, + 0.249286745170910, + 0.873821971016996, + 0.063089014491502, + 0.063089014491502, + 0.053145049844817, + 0.310352451033784, + 0.636502499121399, + 0.310352451033784, + 0.636502499121399, + 0.053145049844817, +] + +const trianglequadDunavant12b = @SVector [ + 0.249286745170910, + 0.249286745170910, + 0.501426509658179, + 0.063089014491502, + 0.063089014491502, + 0.873821971016996, + 0.310352451033784, + 0.636502499121399, + 0.053145049844817, + 0.053145049844817, + 0.310352451033784, + 0.636502499121399, +] + +const trianglequadDunavant12w = @SVector [ + 0.116786275726379, + 0.116786275726379, + 0.116786275726379, + 0.050844906370207, + 0.050844906370207, + 0.050844906370207, + 0.082851075618374, + 0.082851075618374, + 0.082851075618374, + 0.082851075618374, + 0.082851075618374, + 0.082851075618374, +] + +const trianglequadDunavant13a = @SVector [ + 0.333333333333333, + 0.479308067841920, + 0.260345966079040, + 0.260345966079040, + 0.869739794195568, + 0.065130102902216, + 0.065130102902216, + 0.048690315425316, + 0.312865496004874, + 0.638444188569810, + 0.312865496004874, + 0.638444188569810, + 0.048690315425316, +] + +const trianglequadDunavant13b = @SVector [ + 0.333333333333333, + 0.260345966079040, + 0.260345966079040, + 0.479308067841920, + 0.065130102902216, + 0.065130102902216, + 0.869739794195568, + 0.312865496004874, + 0.638444188569810, + 0.048690315425316, + 0.048690315425316, + 0.312865496004874, + 0.638444188569810, +] + +const trianglequadDunavant13w = @SVector [ + -0.149570044467682, + 0.175615257433208, + 0.175615257433208, + 0.175615257433208, + 0.053347235608838, + 0.053347235608838, + 0.053347235608838, + 0.077113760890257, + 0.077113760890257, + 0.077113760890257, + 0.077113760890257, + 0.077113760890257, + 0.077113760890257, +] + +const trianglequadDunavant16a = @SVector [ + 0.333333333333333, + 0.081414823414554, + 0.459292588292723, + 0.459292588292723, + 0.658861384496480, + 0.170569307751760, + 0.170569307751760, + 0.898905543365938, + 0.050547228317031, + 0.050547228317031, + 0.008394777409958, + 0.263112829634638, + 0.728492392955404, + 0.263112829634638, + 0.728492392955404, + 0.008394777409958, +] + +const trianglequadDunavant16b = @SVector [ + 0.333333333333333, + 0.459292588292723, + 0.459292588292723, + 0.081414823414554, + 0.170569307751760, + 0.170569307751760, + 0.658861384496480, + 0.050547228317031, + 0.050547228317031, + 0.898905543365938, + 0.263112829634638, + 0.728492392955404, + 0.008394777409958, + 0.008394777409958, + 0.263112829634638, + 0.728492392955404, +] + +const trianglequadDunavant16w = @SVector [ + 0.144315607677787, + 0.095091634267285, + 0.095091634267285, + 0.095091634267285, + 0.103217370534718, + 0.103217370534718, + 0.103217370534718, + 0.032458497623198, + 0.032458497623198, + 0.032458497623198, + 0.027230314174435, + 0.027230314174435, + 0.027230314174435, + 0.027230314174435, + 0.027230314174435, + 0.027230314174435, +] + +const trianglequadDunavant19a = @SVector [ + 0.333333333333333, + 0.020634961602525, + 0.489682519198738, + 0.489682519198738, + 0.125820817014127, + 0.437089591492937, + 0.437089591492937, + 0.623592928761935, + 0.188203535619033, + 0.188203535619033, + 0.910540973211095, + 0.044729513394453, + 0.044729513394453, + 0.036838412054736, + 0.221962989160766, + 0.741198598784498, + 0.221962989160766, + 0.741198598784498, + 0.036838412054736, +] + +const trianglequadDunavant19b = @SVector [ + 0.333333333333333, + 0.489682519198738, + 0.489682519198738, + 0.020634961602525, + 0.437089591492937, + 0.437089591492937, + 0.125820817014127, + 0.188203535619033, + 0.188203535619033, + 0.623592928761935, + 0.044729513394453, + 0.044729513394453, + 0.910540973211095, + 0.221962989160766, + 0.741198598784498, + 0.036838412054736, + 0.036838412054736, + 0.221962989160766, + 0.741198598784498, +] + +const trianglequadDunavant19w = @SVector [ + 0.097135796282799, + 0.031334700227139, + 0.031334700227139, + 0.031334700227139, + 0.077827541004774, + 0.077827541004774, + 0.077827541004774, + 0.079647738927210, + 0.079647738927210, + 0.079647738927210, + 0.025577675658698, + 0.025577675658698, + 0.025577675658698, + 0.043283539377289, + 0.043283539377289, + 0.043283539377289, + 0.043283539377289, + 0.043283539377289, + 0.043283539377289, +] + +const trianglequadDunavant25a = @SVector [ + 0.333333333333333, + 0.028844733232685, + 0.485577633383657, + 0.485577633383657, + 0.781036849029926, + 0.109481575485037, + 0.109481575485037, + 0.141707219414880, + 0.307939838764121, + 0.550352941820999, + 0.307939838764121, + 0.550352941820999, + 0.141707219414880, + 0.025003534762686, + 0.246672560639903, + 0.728323904597411, + 0.246672560639903, + 0.728323904597411, + 0.025003534762686, + 0.009540815400299, + 0.066803251012200, + 0.923655933587500, + 0.066803251012200, + 0.923655933587500, + 0.009540815400299, +] + +const trianglequadDunavant25b = @SVector [ + 0.333333333333333, + 0.485577633383657, + 0.485577633383657, + 0.028844733232685, + 0.109481575485037, + 0.109481575485037, + 0.781036849029926, + 0.307939838764121, + 0.550352941820999, + 0.141707219414880, + 0.141707219414880, + 0.307939838764121, + 0.550352941820999, + 0.246672560639903, + 0.728323904597411, + 0.025003534762686, + 0.025003534762686, + 0.246672560639903, + 0.728323904597411, + 0.066803251012200, + 0.923655933587500, + 0.009540815400299, + 0.009540815400299, + 0.066803251012200, + 0.923655933587500, +] + +const trianglequadDunavant25w = @SVector [ + 0.090817990382754, + 0.036725957756467, + 0.036725957756467, + 0.036725957756467, + 0.045321059435528, + 0.045321059435528, + 0.045321059435528, + 0.072757916845420, + 0.072757916845420, + 0.072757916845420, + 0.072757916845420, + 0.072757916845420, + 0.072757916845420, + 0.028327242531057, + 0.028327242531057, + 0.028327242531057, + 0.028327242531057, + 0.028327242531057, + 0.028327242531057, + 0.009421666963733, + 0.009421666963733, + 0.009421666963733, + 0.009421666963733, + 0.009421666963733, + 0.009421666963733, +] + +const trianglequadDunavant27a = @SVector [ + -0.069222096541517, + 0.534611048270758, + 0.534611048270758, + 0.202061394068290, + 0.398969302965855, + 0.398969302965855, + 0.593380199137435, + 0.203309900431282, + 0.203309900431282, + 0.761298175434837, + 0.119350912282581, + 0.119350912282581, + 0.935270103777448, + 0.032364948111276, + 0.032364948111276, + 0.050178138310495, + 0.356620648261293, + 0.593201213428213, + 0.356620648261293, + 0.593201213428213, + 0.050178138310495, + 0.021022016536166, + 0.171488980304042, + 0.807489003159792, + 0.171488980304042, + 0.807489003159792, + 0.021022016536166, +] + +const trianglequadDunavant27b = @SVector [ + 0.534611048270758, + 0.534611048270758, + -0.069222096541517, + 0.398969302965855, + 0.398969302965855, + 0.202061394068290, + 0.203309900431282, + 0.203309900431282, + 0.593380199137435, + 0.119350912282581, + 0.119350912282581, + 0.761298175434837, + 0.032364948111276, + 0.032364948111276, + 0.935270103777448, + 0.356620648261293, + 0.593201213428213, + 0.050178138310495, + 0.050178138310495, + 0.356620648261293, + 0.593201213428213, + 0.171488980304042, + 0.807489003159792, + 0.021022016536166, + 0.021022016536166, + 0.171488980304042, + 0.807489003159792, +] + +const trianglequadDunavant27w = @SVector [ + 0.000927006328961, + 0.000927006328961, + 0.000927006328961, + 0.077149534914813, + 0.077149534914813, + 0.077149534914813, + 0.059322977380774, + 0.059322977380774, + 0.059322977380774, + 0.036184540503418, + 0.036184540503418, + 0.036184540503418, + 0.013659731002678, + 0.013659731002678, + 0.013659731002678, + 0.052337111962204, + 0.052337111962204, + 0.052337111962204, + 0.052337111962204, + 0.052337111962204, + 0.052337111962204, + 0.020707659639141, + 0.020707659639141, + 0.020707659639141, + 0.020707659639141, + 0.020707659639141, + 0.020707659639141, +] + +const trianglequadDunavant33a = @SVector [ + 0.023565220452390, + 0.488217389773805, + 0.488217389773805, + 0.120551215411079, + 0.439724392294460, + 0.439724392294460, + 0.457579229975768, + 0.271210385012116, + 0.271210385012116, + 0.744847708916828, + 0.127576145541586, + 0.127576145541586, + 0.957365299093579, + 0.021317350453210, + 0.021317350453210, + 0.115343494534698, + 0.275713269685514, + 0.608943235779788, + 0.275713269685514, + 0.608943235779788, + 0.115343494534698, + 0.022838332222257, + 0.281325580989940, + 0.695836086787803, + 0.281325580989940, + 0.695836086787803, + 0.022838332222257, + 0.025734050548330, + 0.116251915907597, + 0.858014033544073, + 0.116251915907597, + 0.858014033544073, + 0.025734050548330, +] + +const trianglequadDunavant33b = @SVector [ + 0.488217389773805, + 0.488217389773805, + 0.023565220452390, + 0.439724392294460, + 0.439724392294460, + 0.120551215411079, + 0.271210385012116, + 0.271210385012116, + 0.457579229975768, + 0.127576145541586, + 0.127576145541586, + 0.744847708916828, + 0.021317350453210, + 0.021317350453210, + 0.957365299093579, + 0.275713269685514, + 0.608943235779788, + 0.115343494534698, + 0.115343494534698, + 0.275713269685514, + 0.608943235779788, + 0.281325580989940, + 0.695836086787803, + 0.022838332222257, + 0.022838332222257, + 0.281325580989940, + 0.695836086787803, + 0.116251915907597, + 0.858014033544073, + 0.025734050548330, + 0.025734050548330, + 0.116251915907597, + 0.858014033544073, +] + +const trianglequadDunavant33w = @SVector [ + 0.025731066440455, + 0.025731066440455, + 0.025731066440455, + 0.043692544538038, + 0.043692544538038, + 0.043692544538038, + 0.062858224217885, + 0.062858224217885, + 0.062858224217885, + 0.034796112930709, + 0.034796112930709, + 0.034796112930709, + 0.006166261051559, + 0.006166261051559, + 0.006166261051559, + 0.040371557766381, + 0.040371557766381, + 0.040371557766381, + 0.040371557766381, + 0.040371557766381, + 0.040371557766381, + 0.022356773202303, + 0.022356773202303, + 0.022356773202303, + 0.022356773202303, + 0.022356773202303, + 0.022356773202303, + 0.017316231108659, + 0.017316231108659, + 0.017316231108659, + 0.017316231108659, + 0.017316231108659, + 0.017316231108659, +] + +const trianglequadDunavant37a = @SVector [ + 0.333333333333333, + 0.009903630120591, + 0.495048184939705, + 0.495048184939705, + 0.062566729780852, + 0.468716635109574, + 0.468716635109574, + 0.170957326397447, + 0.414521336801277, + 0.414521336801277, + 0.541200855914337, + 0.229399572042831, + 0.229399572042831, + 0.771151009607340, + 0.114424495196330, + 0.114424495196330, + 0.950377217273082, + 0.024811391363459, + 0.024811391363459, + 0.094853828379579, + 0.268794997058761, + 0.636351174561660, + 0.268794997058761, + 0.636351174561660, + 0.094853828379579, + 0.018100773278807, + 0.291730066734288, + 0.690169159986905, + 0.291730066734288, + 0.690169159986905, + 0.018100773278807, + 0.022233076674090, + 0.126357385491669, + 0.851409537834241, + 0.126357385491669, + 0.851409537834241, + 0.022233076674090, +] + +const trianglequadDunavant37b = @SVector [ + 0.333333333333333, + 0.495048184939705, + 0.495048184939705, + 0.009903630120591, + 0.468716635109574, + 0.468716635109574, + 0.062566729780852, + 0.414521336801277, + 0.414521336801277, + 0.170957326397447, + 0.229399572042831, + 0.229399572042831, + 0.541200855914337, + 0.114424495196330, + 0.114424495196330, + 0.771151009607340, + 0.024811391363459, + 0.024811391363459, + 0.950377217273082, + 0.268794997058761, + 0.636351174561660, + 0.094853828379579, + 0.094853828379579, + 0.268794997058761, + 0.636351174561660, + 0.291730066734288, + 0.690169159986905, + 0.018100773278807, + 0.018100773278807, + 0.291730066734288, + 0.690169159986905, + 0.126357385491669, + 0.851409537834241, + 0.022233076674090, + 0.022233076674090, + 0.126357385491669, + 0.851409537834241, +] + +const trianglequadDunavant37w = @SVector [ + 0.052520923400802, + 0.011280145209330, + 0.011280145209330, + 0.011280145209330, + 0.031423518362454, + 0.031423518362454, + 0.031423518362454, + 0.047072502504194, + 0.047072502504194, + 0.047072502504194, + 0.047363586536355, + 0.047363586536355, + 0.047363586536355, + 0.031167529045794, + 0.031167529045794, + 0.031167529045794, + 0.007975771465074, + 0.007975771465074, + 0.007975771465074, + 0.036848402728732, + 0.036848402728732, + 0.036848402728732, + 0.036848402728732, + 0.036848402728732, + 0.036848402728732, + 0.017401463303822, + 0.017401463303822, + 0.017401463303822, + 0.017401463303822, + 0.017401463303822, + 0.017401463303822, + 0.015521786839045, + 0.015521786839045, + 0.015521786839045, + 0.015521786839045, + 0.015521786839045, + 0.015521786839045, +] + +const trianglequadDunavant42a = @SVector [ + 0.022072179275643, + 0.488963910362179, + 0.488963910362179, + 0.164710561319092, + 0.417644719340454, + 0.417644719340454, + 0.453044943382323, + 0.273477528308839, + 0.273477528308839, + 0.645588935174913, + 0.177205532412543, + 0.177205532412543, + 0.876400233818255, + 0.061799883090873, + 0.061799883090873, + 0.961218077502598, + 0.019390961248701, + 0.019390961248701, + 0.057124757403648, + 0.172266687821356, + 0.770608554774996, + 0.172266687821356, + 0.770608554774996, + 0.057124757403648, + 0.092916249356972, + 0.336861459796345, + 0.570222290846683, + 0.336861459796345, + 0.570222290846683, + 0.092916249356972, + 0.014646950055654, + 0.298372882136258, + 0.686980167808088, + 0.298372882136258, + 0.686980167808088, + 0.014646950055654, + 0.001268330932872, + 0.118974497696957, + 0.879757171370171, + 0.118974497696957, + 0.879757171370171, + 0.001268330932872, +] + +const trianglequadDunavant42b = @SVector [ + 0.488963910362179, + 0.488963910362179, + 0.022072179275643, + 0.417644719340454, + 0.417644719340454, + 0.164710561319092, + 0.273477528308839, + 0.273477528308839, + 0.453044943382323, + 0.177205532412543, + 0.177205532412543, + 0.645588935174913, + 0.061799883090873, + 0.061799883090873, + 0.876400233818255, + 0.019390961248701, + 0.019390961248701, + 0.961218077502598, + 0.172266687821356, + 0.770608554774996, + 0.057124757403648, + 0.057124757403648, + 0.172266687821356, + 0.770608554774996, + 0.336861459796345, + 0.570222290846683, + 0.092916249356972, + 0.092916249356972, + 0.336861459796345, + 0.570222290846683, + 0.298372882136258, + 0.686980167808088, + 0.014646950055654, + 0.014646950055654, + 0.298372882136258, + 0.686980167808088, + 0.118974497696957, + 0.879757171370171, + 0.001268330932872, + 0.001268330932872, + 0.118974497696957, + 0.879757171370171, +] + +const trianglequadDunavant42w = @SVector [ + 0.021883581369429, + 0.021883581369429, + 0.021883581369429, + 0.032788353544125, + 0.032788353544125, + 0.032788353544125, + 0.051774104507292, + 0.051774104507292, + 0.051774104507292, + 0.042162588736993, + 0.042162588736993, + 0.042162588736993, + 0.014433699669777, + 0.014433699669777, + 0.014433699669777, + 0.004923403602400, + 0.004923403602400, + 0.004923403602400, + 0.024665753212564, + 0.024665753212564, + 0.024665753212564, + 0.024665753212564, + 0.024665753212564, + 0.024665753212564, + 0.038571510787061, + 0.038571510787061, + 0.038571510787061, + 0.038571510787061, + 0.038571510787061, + 0.038571510787061, + 0.014436308113534, + 0.014436308113534, + 0.014436308113534, + 0.014436308113534, + 0.014436308113534, + 0.014436308113534, + 0.005010228838501, + 0.005010228838501, + 0.005010228838501, + 0.005010228838501, + 0.005010228838501, + 0.005010228838501, +] + +const trianglequadDunavant48a = @SVector [ + -0.013945833716486, + 0.506972916858243, + 0.506972916858243, + 0.137187291433955, + 0.431406354283023, + 0.431406354283023, + 0.444612710305711, + 0.277693644847144, + 0.277693644847144, + 0.747070217917492, + 0.126464891041254, + 0.126464891041254, + 0.858383228050628, + 0.070808385974686, + 0.070808385974686, + 0.962069659517853, + 0.018965170241073, + 0.018965170241073, + 0.133734161966621, + 0.261311371140087, + 0.604954466893291, + 0.261311371140087, + 0.604954466893291, + 0.133734161966621, + 0.036366677396917, + 0.388046767090269, + 0.575586555512814, + 0.388046767090269, + 0.575586555512814, + 0.036366677396917, + -0.010174883126571, + 0.285712220049916, + 0.724462663076655, + 0.285712220049916, + 0.724462663076655, + -0.010174883126571, + 0.036843869875878, + 0.215599664072284, + 0.747556466051838, + 0.215599664072284, + 0.747556466051838, + 0.036843869875878, + 0.012459809331199, + 0.103575616576386, + 0.883964574092416, + 0.103575616576386, + 0.883964574092416, + 0.012459809331199, +] + +const trianglequadDunavant48b = @SVector [ + 0.506972916858243, + 0.506972916858243, + -0.013945833716486, + 0.431406354283023, + 0.431406354283023, + 0.137187291433955, + 0.277693644847144, + 0.277693644847144, + 0.444612710305711, + 0.126464891041254, + 0.126464891041254, + 0.747070217917492, + 0.070808385974686, + 0.070808385974686, + 0.858383228050628, + 0.018965170241073, + 0.018965170241073, + 0.962069659517853, + 0.261311371140087, + 0.604954466893291, + 0.133734161966621, + 0.133734161966621, + 0.261311371140087, + 0.604954466893291, + 0.388046767090269, + 0.575586555512814, + 0.036366677396917, + 0.036366677396917, + 0.388046767090269, + 0.575586555512814, + 0.285712220049916, + 0.724462663076655, + -0.010174883126571, + -0.010174883126571, + 0.285712220049916, + 0.724462663076655, + 0.215599664072284, + 0.747556466051838, + 0.036843869875878, + 0.036843869875878, + 0.215599664072284, + 0.747556466051838, + 0.103575616576386, + 0.883964574092416, + 0.012459809331199, + 0.012459809331199, + 0.103575616576386, + 0.883964574092416, +] + +const trianglequadDunavant48w = @SVector [ + 0.001916875642849, + 0.001916875642849, + 0.001916875642849, + 0.044249027271145, + 0.044249027271145, + 0.044249027271145, + 0.051186548718852, + 0.051186548718852, + 0.051186548718852, + 0.023687735870688, + 0.023687735870688, + 0.023687735870688, + 0.013289775690021, + 0.013289775690021, + 0.013289775690021, + 0.004748916608192, + 0.004748916608192, + 0.004748916608192, + 0.038550072599593, + 0.038550072599593, + 0.038550072599593, + 0.038550072599593, + 0.038550072599593, + 0.038550072599593, + 0.027215814320624, + 0.027215814320624, + 0.027215814320624, + 0.027215814320624, + 0.027215814320624, + 0.027215814320624, + 0.002182077366797, + 0.002182077366797, + 0.002182077366797, + 0.002182077366797, + 0.002182077366797, + 0.002182077366797, + 0.021505319847731, + 0.021505319847731, + 0.021505319847731, + 0.021505319847731, + 0.021505319847731, + 0.021505319847731, + 0.007673942631049, + 0.007673942631049, + 0.007673942631049, + 0.007673942631049, + 0.007673942631049, + 0.007673942631049, +] + +const trianglequadDunavant52a = @SVector [ + 0.333333333333333, + 0.005238916103123, + 0.497380541948438, + 0.497380541948438, + 0.173061122901295, + 0.413469438549352, + 0.413469438549352, + 0.059082801866017, + 0.470458599066991, + 0.470458599066991, + 0.518892500060958, + 0.240553749969521, + 0.240553749969521, + 0.704068411554854, + 0.147965794222573, + 0.147965794222573, + 0.849069624685052, + 0.075465187657474, + 0.075465187657474, + 0.966807194753950, + 0.016596402623025, + 0.016596402623025, + 0.103575692245252, + 0.296555596579887, + 0.599868711174861, + 0.296555596579887, + 0.599868711174861, + 0.103575692245252, + 0.020083411655416, + 0.337723063403079, + 0.642193524941505, + 0.337723063403079, + 0.642193524941505, + 0.020083411655416, + -0.004341002614139, + 0.204748281642812, + 0.799592720971327, + 0.204748281642812, + 0.799592720971327, + -0.004341002614139, + 0.041941786468010, + 0.189358492130623, + 0.768699721401368, + 0.189358492130623, + 0.768699721401368, + 0.041941786468010, + 0.014317320230681, + 0.085283615682657, + 0.900399064086661, + 0.085283615682657, + 0.900399064086661, + 0.014317320230681, +] + +const trianglequadDunavant52b = @SVector [ + 0.333333333333333, + 0.497380541948438, + 0.497380541948438, + 0.005238916103123, + 0.413469438549352, + 0.413469438549352, + 0.173061122901295, + 0.470458599066991, + 0.470458599066991, + 0.059082801866017, + 0.240553749969521, + 0.240553749969521, + 0.518892500060958, + 0.147965794222573, + 0.147965794222573, + 0.704068411554854, + 0.075465187657474, + 0.075465187657474, + 0.849069624685052, + 0.016596402623025, + 0.016596402623025, + 0.966807194753950, + 0.296555596579887, + 0.599868711174861, + 0.103575692245252, + 0.103575692245252, + 0.296555596579887, + 0.599868711174861, + 0.337723063403079, + 0.642193524941505, + 0.020083411655416, + 0.020083411655416, + 0.337723063403079, + 0.642193524941505, + 0.204748281642812, + 0.799592720971327, + -0.004341002614139, + -0.004341002614139, + 0.204748281642812, + 0.799592720971327, + 0.189358492130623, + 0.768699721401368, + 0.041941786468010, + 0.041941786468010, + 0.189358492130623, + 0.768699721401368, + 0.085283615682657, + 0.900399064086661, + 0.014317320230681, + 0.014317320230681, + 0.085283615682657, + 0.900399064086661, +] + +const trianglequadDunavant52w = @SVector [ + 0.046875697427642, + 0.006405878578585, + 0.006405878578585, + 0.006405878578585, + 0.041710296739387, + 0.041710296739387, + 0.041710296739387, + 0.026891484250064, + 0.026891484250064, + 0.026891484250064, + 0.042132522761650, + 0.042132522761650, + 0.042132522761650, + 0.030000266842773, + 0.030000266842773, + 0.030000266842773, + 0.014200098925024, + 0.014200098925024, + 0.014200098925024, + 0.003582462351273, + 0.003582462351273, + 0.003582462351273, + 0.032773147460627, + 0.032773147460627, + 0.032773147460627, + 0.032773147460627, + 0.032773147460627, + 0.032773147460627, + 0.015298306248441, + 0.015298306248441, + 0.015298306248441, + 0.015298306248441, + 0.015298306248441, + 0.015298306248441, + 0.002386244192839, + 0.002386244192839, + 0.002386244192839, + 0.002386244192839, + 0.002386244192839, + 0.002386244192839, + 0.019084792755899, + 0.019084792755899, + 0.019084792755899, + 0.019084792755899, + 0.019084792755899, + 0.019084792755899, + 0.006850054546542, + 0.006850054546542, + 0.006850054546542, + 0.006850054546542, + 0.006850054546542, + 0.006850054546542, +] + +const trianglequadDunavant61a = @SVector [ + 0.333333333333333, + 0.005658918886452, + 0.497170540556774, + 0.497170540556774, + 0.035647354750751, + 0.482176322624625, + 0.482176322624625, + 0.099520061958437, + 0.450239969020782, + 0.450239969020782, + 0.199467521245206, + 0.400266239377397, + 0.400266239377397, + 0.495717464058095, + 0.252141267970953, + 0.252141267970953, + 0.675905990683077, + 0.162047004658461, + 0.162047004658461, + 0.848248235478508, + 0.075875882260746, + 0.075875882260746, + 0.968690546064356, + 0.015654726967822, + 0.015654726967822, + 0.010186928826919, + 0.334319867363658, + 0.655493203809423, + 0.334319867363658, + 0.655493203809423, + 0.010186928826919, + 0.135440871671036, + 0.292221537796944, + 0.572337590532020, + 0.292221537796944, + 0.572337590532020, + 0.135440871671036, + 0.054423924290583, + 0.319574885423190, + 0.626001190286228, + 0.319574885423190, + 0.626001190286228, + 0.054423924290583, + 0.012868560833637, + 0.190704224192292, + 0.796427214974071, + 0.190704224192292, + 0.796427214974071, + 0.012868560833637, + 0.067165782413524, + 0.180483211648746, + 0.752351005937729, + 0.180483211648746, + 0.752351005937729, + 0.067165782413524, + 0.014663182224828, + 0.080711313679564, + 0.904625504095608, + 0.080711313679564, + 0.904625504095608, + 0.014663182224828, +] + +const trianglequadDunavant61b = @SVector [ + 0.333333333333333, + 0.497170540556774, + 0.497170540556774, + 0.005658918886452, + 0.482176322624625, + 0.482176322624625, + 0.035647354750751, + 0.450239969020782, + 0.450239969020782, + 0.099520061958437, + 0.400266239377397, + 0.400266239377397, + 0.199467521245206, + 0.252141267970953, + 0.252141267970953, + 0.495717464058095, + 0.162047004658461, + 0.162047004658461, + 0.675905990683077, + 0.075875882260746, + 0.075875882260746, + 0.848248235478508, + 0.015654726967822, + 0.015654726967822, + 0.968690546064356, + 0.334319867363658, + 0.655493203809423, + 0.010186928826919, + 0.010186928826919, + 0.334319867363658, + 0.655493203809423, + 0.292221537796944, + 0.572337590532020, + 0.135440871671036, + 0.135440871671036, + 0.292221537796944, + 0.572337590532020, + 0.319574885423190, + 0.626001190286228, + 0.054423924290583, + 0.054423924290583, + 0.319574885423190, + 0.626001190286228, + 0.190704224192292, + 0.796427214974071, + 0.012868560833637, + 0.012868560833637, + 0.190704224192292, + 0.796427214974071, + 0.180483211648746, + 0.752351005937729, + 0.067165782413524, + 0.067165782413524, + 0.180483211648746, + 0.752351005937729, + 0.080711313679564, + 0.904625504095608, + 0.014663182224828, + 0.014663182224828, + 0.080711313679564, + 0.904625504095608, +] + +const trianglequadDunavant61w = @SVector [ + 0.033437199290803, + 0.005093415440507, + 0.005093415440507, + 0.005093415440507, + 0.014670864527638, + 0.014670864527638, + 0.014670864527638, + 0.024350878353672, + 0.024350878353672, + 0.024350878353672, + 0.031107550868969, + 0.031107550868969, + 0.031107550868969, + 0.031257111218620, + 0.031257111218620, + 0.031257111218620, + 0.024815654339665, + 0.024815654339665, + 0.024815654339665, + 0.014056073070557, + 0.014056073070557, + 0.014056073070557, + 0.003194676173779, + 0.003194676173779, + 0.003194676173779, + 0.008119655318993, + 0.008119655318993, + 0.008119655318993, + 0.008119655318993, + 0.008119655318993, + 0.008119655318993, + 0.026805742283163, + 0.026805742283163, + 0.026805742283163, + 0.026805742283163, + 0.026805742283163, + 0.026805742283163, + 0.018459993210822, + 0.018459993210822, + 0.018459993210822, + 0.018459993210822, + 0.018459993210822, + 0.018459993210822, + 0.008476868534328, + 0.008476868534328, + 0.008476868534328, + 0.008476868534328, + 0.008476868534328, + 0.008476868534328, + 0.018292796770025, + 0.018292796770025, + 0.018292796770025, + 0.018292796770025, + 0.018292796770025, + 0.018292796770025, + 0.006665632004165, + 0.006665632004165, + 0.006665632004165, + 0.006665632004165, + 0.006665632004165, + 0.006665632004165, +] + +const trianglequadDunavant70a = @SVector [ + 0.333333333333333, + 0.013310382738157, + 0.493344808630921, + 0.493344808630921, + 0.061578811516086, + 0.469210594241957, + 0.469210594241957, + 0.127437208225989, + 0.436281395887006, + 0.436281395887006, + 0.210307658653168, + 0.394846170673416, + 0.394846170673416, + 0.500410862393686, + 0.249794568803157, + 0.249794568803157, + 0.677135612512315, + 0.161432193743843, + 0.161432193743843, + 0.846803545029257, + 0.076598227485371, + 0.076598227485371, + 0.951495121293100, + 0.024252439353450, + 0.024252439353450, + 0.913707265566071, + 0.043146367216965, + 0.043146367216965, + 0.008430536202420, + 0.358911494940944, + 0.632657968856636, + 0.358911494940944, + 0.632657968856636, + 0.008430536202420, + 0.131186551737188, + 0.294402476751957, + 0.574410971510855, + 0.294402476751957, + 0.574410971510855, + 0.131186551737188, + 0.050203151565675, + 0.325017801641814, + 0.624779046792512, + 0.325017801641814, + 0.624779046792512, + 0.050203151565675, + 0.066329263810916, + 0.184737559666046, + 0.748933176523037, + 0.184737559666046, + 0.748933176523037, + 0.066329263810916, + 0.011996194566236, + 0.218796800013321, + 0.769207005420443, + 0.218796800013321, + 0.769207005420443, + 0.011996194566236, + 0.014858100590125, + 0.101179597136408, + 0.883962302273467, + 0.101179597136408, + 0.883962302273467, + 0.014858100590125, + -0.035222015287949, + 0.020874755282586, + 1.014347260005363, + 0.020874755282586, + 1.014347260005363, + -0.035222015287949, +] + +const trianglequadDunavant70b = @SVector [ + 0.333333333333333, + 0.493344808630921, + 0.493344808630921, + 0.013310382738157, + 0.469210594241957, + 0.469210594241957, + 0.061578811516086, + 0.436281395887006, + 0.436281395887006, + 0.127437208225989, + 0.394846170673416, + 0.394846170673416, + 0.210307658653168, + 0.249794568803157, + 0.249794568803157, + 0.500410862393686, + 0.161432193743843, + 0.161432193743843, + 0.677135612512315, + 0.076598227485371, + 0.076598227485371, + 0.846803545029257, + 0.024252439353450, + 0.024252439353450, + 0.951495121293100, + 0.043146367216965, + 0.043146367216965, + 0.913707265566071, + 0.358911494940944, + 0.632657968856636, + 0.008430536202420, + 0.008430536202420, + 0.358911494940944, + 0.632657968856636, + 0.294402476751957, + 0.574410971510855, + 0.131186551737188, + 0.131186551737188, + 0.294402476751957, + 0.574410971510855, + 0.325017801641814, + 0.624779046792512, + 0.050203151565675, + 0.050203151565675, + 0.325017801641814, + 0.624779046792512, + 0.184737559666046, + 0.748933176523037, + 0.066329263810916, + 0.066329263810916, + 0.184737559666046, + 0.748933176523037, + 0.218796800013321, + 0.769207005420443, + 0.011996194566236, + 0.011996194566236, + 0.218796800013321, + 0.769207005420443, + 0.101179597136408, + 0.883962302273467, + 0.014858100590125, + 0.014858100590125, + 0.101179597136408, + 0.883962302273467, + 0.020874755282586, + 1.014347260005363, + -0.035222015287949, + -0.035222015287949, + 0.020874755282586, + 1.014347260005363, +] + +const trianglequadDunavant70w = @SVector [ + 0.030809939937647, + 0.009072436679404, + 0.009072436679404, + 0.009072436679404, + 0.018761316939594, + 0.018761316939594, + 0.018761316939594, + 0.019441097985477, + 0.019441097985477, + 0.019441097985477, + 0.027753948610810, + 0.027753948610810, + 0.027753948610810, + 0.032256225351457, + 0.032256225351457, + 0.032256225351457, + 0.025074032616922, + 0.025074032616922, + 0.025074032616922, + 0.015271927971832, + 0.015271927971832, + 0.015271927971832, + 0.006793922022963, + 0.006793922022963, + 0.006793922022963, + -0.002223098729920, + -0.002223098729920, + -0.002223098729920, + 0.006331914076406, + 0.006331914076406, + 0.006331914076406, + 0.006331914076406, + 0.006331914076406, + 0.006331914076406, + 0.027257538049138, + 0.027257538049138, + 0.027257538049138, + 0.027257538049138, + 0.027257538049138, + 0.027257538049138, + 0.017676785649465, + 0.017676785649465, + 0.017676785649465, + 0.017676785649465, + 0.017676785649465, + 0.017676785649465, + 0.018379484638070, + 0.018379484638070, + 0.018379484638070, + 0.018379484638070, + 0.018379484638070, + 0.018379484638070, + 0.008104732808192, + 0.008104732808192, + 0.008104732808192, + 0.008104732808192, + 0.008104732808192, + 0.008104732808192, + 0.007634129070725, + 0.007634129070725, + 0.007634129070725, + 0.007634129070725, + 0.007634129070725, + 0.007634129070725, + 0.000046187660794, + 0.000046187660794, + 0.000046187660794, + 0.000046187660794, + 0.000046187660794, + 0.000046187660794, +] + +const trianglequadDunavant73a = @SVector [ + 0.333333333333333, + 0.020780025853987, + 0.489609987073006, + 0.489609987073006, + 0.090926214604215, + 0.454536892697893, + 0.454536892697893, + 0.197166638701138, + 0.401416680649431, + 0.401416680649431, + 0.488896691193805, + 0.255551654403098, + 0.255551654403098, + 0.645844115695741, + 0.177077942152130, + 0.177077942152130, + 0.779877893544096, + 0.110061053227952, + 0.110061053227952, + 0.888942751496321, + 0.055528624251840, + 0.055528624251840, + 0.974756272445543, + 0.012621863777229, + 0.012621863777229, + 0.003611417848412, + 0.395754787356943, + 0.600633794794645, + 0.395754787356943, + 0.600633794794645, + 0.003611417848412, + 0.134466754530780, + 0.307929983880436, + 0.557603261588784, + 0.307929983880436, + 0.557603261588784, + 0.134466754530780, + 0.014446025776115, + 0.264566948406520, + 0.720987025817365, + 0.264566948406520, + 0.720987025817365, + 0.014446025776115, + 0.046933578838178, + 0.358539352205951, + 0.594527068955871, + 0.358539352205951, + 0.594527068955871, + 0.046933578838178, + 0.002861120350567, + 0.157807405968595, + 0.839331473680839, + 0.157807405968595, + 0.839331473680839, + 0.002861120350567, + 0.223861424097916, + 0.075050596975911, + 0.701087978926173, + 0.075050596975911, + 0.701087978926173, + 0.223861424097916, + 0.034647074816760, + 0.142421601113383, + 0.822931324069857, + 0.142421601113383, + 0.822931324069857, + 0.034647074816760, + 0.010161119296278, + 0.065494628082938, + 0.924344252620784, + 0.065494628082938, + 0.924344252620784, + 0.010161119296278, +] + +const trianglequadDunavant73b = @SVector [ + 0.333333333333333, + 0.489609987073006, + 0.489609987073006, + 0.020780025853987, + 0.454536892697893, + 0.454536892697893, + 0.090926214604215, + 0.401416680649431, + 0.401416680649431, + 0.197166638701138, + 0.255551654403098, + 0.255551654403098, + 0.488896691193805, + 0.177077942152130, + 0.177077942152130, + 0.645844115695741, + 0.110061053227952, + 0.110061053227952, + 0.779877893544096, + 0.055528624251840, + 0.055528624251840, + 0.888942751496321, + 0.012621863777229, + 0.012621863777229, + 0.974756272445543, + 0.395754787356943, + 0.600633794794645, + 0.003611417848412, + 0.003611417848412, + 0.395754787356943, + 0.600633794794645, + 0.307929983880436, + 0.557603261588784, + 0.134466754530780, + 0.134466754530780, + 0.307929983880436, + 0.557603261588784, + 0.264566948406520, + 0.720987025817365, + 0.014446025776115, + 0.014446025776115, + 0.264566948406520, + 0.720987025817365, + 0.358539352205951, + 0.594527068955871, + 0.046933578838178, + 0.046933578838178, + 0.358539352205951, + 0.594527068955871, + 0.157807405968595, + 0.839331473680839, + 0.002861120350567, + 0.002861120350567, + 0.157807405968595, + 0.839331473680839, + 0.075050596975911, + 0.701087978926173, + 0.223861424097916, + 0.223861424097916, + 0.075050596975911, + 0.701087978926173, + 0.142421601113383, + 0.822931324069857, + 0.034647074816760, + 0.034647074816760, + 0.142421601113383, + 0.822931324069857, + 0.065494628082938, + 0.924344252620784, + 0.010161119296278, + 0.010161119296278, + 0.065494628082938, + 0.924344252620784, +] + +const trianglequadDunavant73w = @SVector [ + 0.032906331388919, + 0.010330731891272, + 0.010330731891272, + 0.010330731891272, + 0.022387247263016, + 0.022387247263016, + 0.022387247263016, + 0.030266125869468, + 0.030266125869468, + 0.030266125869468, + 0.030490967802198, + 0.030490967802198, + 0.030490967802198, + 0.024159212741641, + 0.024159212741641, + 0.024159212741641, + 0.016050803586801, + 0.016050803586801, + 0.016050803586801, + 0.008084580261784, + 0.008084580261784, + 0.008084580261784, + 0.002079362027485, + 0.002079362027485, + 0.002079362027485, + 0.003884876904981, + 0.003884876904981, + 0.003884876904981, + 0.003884876904981, + 0.003884876904981, + 0.003884876904981, + 0.025574160612022, + 0.025574160612022, + 0.025574160612022, + 0.025574160612022, + 0.025574160612022, + 0.025574160612022, + 0.008880903573338, + 0.008880903573338, + 0.008880903573338, + 0.008880903573338, + 0.008880903573338, + 0.008880903573338, + 0.016124546761731, + 0.016124546761731, + 0.016124546761731, + 0.016124546761731, + 0.016124546761731, + 0.016124546761731, + 0.002491941817491, + 0.002491941817491, + 0.002491941817491, + 0.002491941817491, + 0.002491941817491, + 0.002491941817491, + 0.018242840118951, + 0.018242840118951, + 0.018242840118951, + 0.018242840118951, + 0.018242840118951, + 0.018242840118951, + 0.010258563736199, + 0.010258563736199, + 0.010258563736199, + 0.010258563736199, + 0.010258563736199, + 0.010258563736199, + 0.003799928855302, + 0.003799928855302, + 0.003799928855302, + 0.003799928855302, + 0.003799928855302, + 0.003799928855302, +] + +const trianglequadDunavant79a = @SVector [ + 0.333333333333333, + -0.001900928704400, + 0.500950464352200, + 0.500950464352200, + 0.023574084130543, + 0.488212957934729, + 0.488212957934729, + 0.089726636099435, + 0.455136681950283, + 0.455136681950283, + 0.196007481363421, + 0.401996259318289, + 0.401996259318289, + 0.488214180481157, + 0.255892909759421, + 0.255892909759421, + 0.647023488009788, + 0.176488255995106, + 0.176488255995106, + 0.791658289326483, + 0.104170855336758, + 0.104170855336758, + 0.893862072318140, + 0.053068963840930, + 0.053068963840930, + 0.916762569607942, + 0.041618715196029, + 0.041618715196029, + 0.976836157186356, + 0.011581921406822, + 0.011581921406822, + 0.048741583664839, + 0.344855770229001, + 0.606402646106160, + 0.344855770229001, + 0.606402646106160, + 0.048741583664839, + 0.006314115948605, + 0.377843269594854, + 0.615842614456541, + 0.377843269594854, + 0.615842614456541, + 0.006314115948605, + 0.134316520547348, + 0.306635479062357, + 0.559048000390295, + 0.306635479062357, + 0.559048000390295, + 0.134316520547348, + 0.013973893962392, + 0.249419362774742, + 0.736606743262866, + 0.249419362774742, + 0.736606743262866, + 0.013973893962392, + 0.075549132909764, + 0.212775724802802, + 0.711675142287434, + 0.212775724802802, + 0.711675142287434, + 0.075549132909764, + -0.008368153208227, + 0.146965436053239, + 0.861402717154987, + 0.146965436053239, + 0.861402717154987, + -0.008368153208227, + 0.026686063258714, + 0.137726978828923, + 0.835586957912363, + 0.137726978828923, + 0.835586957912363, + 0.026686063258714, + 0.010547719294141, + 0.059696109149007, + 0.929756171556853, + 0.059696109149007, + 0.929756171556853, + 0.010547719294141, +] + +const trianglequadDunavant79b = @SVector [ + 0.333333333333333, + 0.500950464352200, + 0.500950464352200, + -0.001900928704400, + 0.488212957934729, + 0.488212957934729, + 0.023574084130543, + 0.455136681950283, + 0.455136681950283, + 0.089726636099435, + 0.401996259318289, + 0.401996259318289, + 0.196007481363421, + 0.255892909759421, + 0.255892909759421, + 0.488214180481157, + 0.176488255995106, + 0.176488255995106, + 0.647023488009788, + 0.104170855336758, + 0.104170855336758, + 0.791658289326483, + 0.053068963840930, + 0.053068963840930, + 0.893862072318140, + 0.041618715196029, + 0.041618715196029, + 0.916762569607942, + 0.011581921406822, + 0.011581921406822, + 0.976836157186356, + 0.344855770229001, + 0.606402646106160, + 0.048741583664839, + 0.048741583664839, + 0.344855770229001, + 0.606402646106160, + 0.377843269594854, + 0.615842614456541, + 0.006314115948605, + 0.006314115948605, + 0.377843269594854, + 0.615842614456541, + 0.306635479062357, + 0.559048000390295, + 0.134316520547348, + 0.134316520547348, + 0.306635479062357, + 0.559048000390295, + 0.249419362774742, + 0.736606743262866, + 0.013973893962392, + 0.013973893962392, + 0.249419362774742, + 0.736606743262866, + 0.212775724802802, + 0.711675142287434, + 0.075549132909764, + 0.075549132909764, + 0.212775724802802, + 0.711675142287434, + 0.146965436053239, + 0.861402717154987, + -0.008368153208227, + -0.008368153208227, + 0.146965436053239, + 0.861402717154987, + 0.137726978828923, + 0.835586957912363, + 0.026686063258714, + 0.026686063258714, + 0.137726978828923, + 0.835586957912363, + 0.059696109149007, + 0.929756171556853, + 0.010547719294141, + 0.010547719294141, + 0.059696109149007, + 0.929756171556853, +] + +const trianglequadDunavant79w = @SVector [ + 0.033057055541624, + 0.000867019185663, + 0.000867019185663, + 0.000867019185663, + 0.011660052716448, + 0.011660052716448, + 0.011660052716448, + 0.022876936356421, + 0.022876936356421, + 0.022876936356421, + 0.030448982673938, + 0.030448982673938, + 0.030448982673938, + 0.030624891725355, + 0.030624891725355, + 0.030624891725355, + 0.024368057676800, + 0.024368057676800, + 0.024368057676800, + 0.015997432032024, + 0.015997432032024, + 0.015997432032024, + 0.007698301815602, + 0.007698301815602, + 0.007698301815602, + -0.000632060497488, + -0.000632060497488, + -0.000632060497488, + 0.001751134301193, + 0.001751134301193, + 0.001751134301193, + 0.016465839189576, + 0.016465839189576, + 0.016465839189576, + 0.016465839189576, + 0.016465839189576, + 0.016465839189576, + 0.004839033540485, + 0.004839033540485, + 0.004839033540485, + 0.004839033540485, + 0.004839033540485, + 0.004839033540485, + 0.025804906534650, + 0.025804906534650, + 0.025804906534650, + 0.025804906534650, + 0.025804906534650, + 0.025804906534650, + 0.008471091054441, + 0.008471091054441, + 0.008471091054441, + 0.008471091054441, + 0.008471091054441, + 0.008471091054441, + 0.018354914106280, + 0.018354914106280, + 0.018354914106280, + 0.018354914106280, + 0.018354914106280, + 0.018354914106280, + 0.000704404677908, + 0.000704404677908, + 0.000704404677908, + 0.000704404677908, + 0.000704404677908, + 0.000704404677908, + 0.010112684927462, + 0.010112684927462, + 0.010112684927462, + 0.010112684927462, + 0.010112684927462, + 0.010112684927462, + 0.003573909385950, + 0.003573909385950, + 0.003573909385950, + 0.003573909385950, + 0.003573909385950, + 0.003573909385950, +] + +const trianglequadDunavantA = @SVector [ + trianglequadDunavant1a, + trianglequadDunavant3a, + trianglequadDunavant4a, + trianglequadDunavant6a, + trianglequadDunavant7a, + trianglequadDunavant12a, + trianglequadDunavant13a, + trianglequadDunavant16a, + trianglequadDunavant19a, + trianglequadDunavant25a, + trianglequadDunavant27a, + trianglequadDunavant33a, + trianglequadDunavant37a, + trianglequadDunavant42a, + trianglequadDunavant48a, + trianglequadDunavant52a, + trianglequadDunavant61a, + trianglequadDunavant70a, + trianglequadDunavant73a, + trianglequadDunavant79a +] + +const trianglequadDunavantB = @SVector [ + trianglequadDunavant1b, + trianglequadDunavant3b, + trianglequadDunavant4b, + trianglequadDunavant6b, + trianglequadDunavant7b, + trianglequadDunavant12b, + trianglequadDunavant13b, + trianglequadDunavant16b, + trianglequadDunavant19b, + trianglequadDunavant25b, + trianglequadDunavant27b, + trianglequadDunavant33b, + trianglequadDunavant37b, + trianglequadDunavant42b, + trianglequadDunavant48b, + trianglequadDunavant52b, + trianglequadDunavant61b, + trianglequadDunavant70b, + trianglequadDunavant73b, + trianglequadDunavant79b +] + +const trianglequadDunavantW = @SVector [ + trianglequadDunavant1w, + trianglequadDunavant3w, + trianglequadDunavant4w, + trianglequadDunavant6w, + trianglequadDunavant7w, + trianglequadDunavant12w, + trianglequadDunavant13w, + trianglequadDunavant16w, + trianglequadDunavant19w, + trianglequadDunavant25w, + trianglequadDunavant27w, + trianglequadDunavant33w, + trianglequadDunavant37w, + trianglequadDunavant42w, + trianglequadDunavant48w, + trianglequadDunavant52w, + trianglequadDunavant61w, + trianglequadDunavant70w, + trianglequadDunavant73w, + trianglequadDunavant79w +] \ No newline at end of file diff --git a/src/quadrature/TriangleGauss.jl b/src/quadrature/triangle/TriangleGauss.jl similarity index 97% rename from src/quadrature/TriangleGauss.jl rename to src/quadrature/triangle/TriangleGauss.jl index 2bbc30d..94bd6cb 100644 --- a/src/quadrature/TriangleGauss.jl +++ b/src/quadrature/triangle/TriangleGauss.jl @@ -1,75 +1,79 @@ -export triangleGaussA, triangleGaussB, triangleGaussW -export trgauss +# First seven cubature rules based on the paper +# Dunavant, D. A. “High Degree Efficient Symmetrical Gaussian Quadrature +# Rules for the Triangle.” International Journal for Numerical Methods +# in Engineering 21, no. 6 (1985): 1129–48. https://doi.org/10.1002/nme.1620210612. -const Trianglegauss1a = [ +# TODO: Add reference to higher order cubature rules + +const trianglequadGauss1a = @SVector [ 0.333333333333333, ] -const Trianglegauss1b = [ +const trianglequadGauss1b = @SVector [ 0.333333333333333, ] -const Trianglegauss1c = [ +const trianglequadGauss1c = @SVector [ 0.333333333333333, ] -const Trianglegauss1w = [ +const trianglequadGauss1w = @SVector [ 1.0, ] -const Trianglegauss3a = [ +const trianglequadGauss3a = @SVector [ 0.666666666666667, 0.166666666666667, 0.166666666666667, ] -const Trianglegauss3b = [ +const trianglequadGauss3b = @SVector [ 0.166666666666667, 0.666666666666667, 0.166666666666667, ] -const Trianglegauss3c = [ +const trianglequadGauss3c = @SVector [ 0.166666666666667, 0.166666666666667, 0.666666666666667, ] -const Trianglegauss3w = [ +const trianglequadGauss3w = @SVector [ 0.333333333333333, 0.333333333333333, 0.333333333333333, ] -const Trianglegauss4a = [ +const trianglequadGauss4a = @SVector [ 0.333333333333333, 0.6, 0.2, 0.2, ] -const Trianglegauss4b = [ +const trianglequadGauss4b = @SVector [ 0.333333333333333, 0.2, 0.6, 0.2, ] -const Trianglegauss4c = [ +const trianglequadGauss4c = @SVector [ 0.333333333333333, 0.2, 0.2, 0.6, ] -const Trianglegauss4w = [ +const trianglequadGauss4w = @SVector [ -0.562500000000000, 0.520833333333333, 0.520833333333333, 0.520833333333333, ] -const Trianglegauss6a = [ +const trianglequadGauss6a = @SVector [ 0.816847572980459, 0.091576213509771, 0.091576213509771, @@ -78,7 +82,7 @@ const Trianglegauss6a = [ 0.445948490915965, ] -const Trianglegauss6b = [ +const trianglequadGauss6b = @SVector [ 0.091576213509771, 0.816847572980459, 0.091576213509771, @@ -87,7 +91,7 @@ const Trianglegauss6b = [ 0.445948490915965, ] -const Trianglegauss6c = [ +const trianglequadGauss6c = @SVector [ 0.091576213509771, 0.091576213509771, 0.816847572980459, @@ -96,7 +100,7 @@ const Trianglegauss6c = [ 0.108103018168070, ] -const Trianglegauss6w = [ +const trianglequadGauss6w = @SVector [ 0.109951743655322, 0.109951743655322, 0.109951743655322, @@ -105,7 +109,7 @@ const Trianglegauss6w = [ 0.223381589678011, ] -const Trianglegauss7a = [ +const trianglequadGauss7a = @SVector [ 0.333333333333333, 0.797426985353087, 0.101286507323456, @@ -115,7 +119,7 @@ const Trianglegauss7a = [ 0.470142064105115, ] -const Trianglegauss7b = [ +const trianglequadGauss7b = @SVector [ 0.333333333333333, 0.101286507323456, 0.797426985353087, @@ -125,7 +129,7 @@ const Trianglegauss7b = [ 0.470142064105115, ] -const Trianglegauss7c = [ +const trianglequadGauss7c = @SVector [ 0.333333333333333, 0.101286507323456, 0.101286507323456, @@ -135,7 +139,7 @@ const Trianglegauss7c = [ 0.059715871789770, ] -const Trianglegauss7w = [ +const trianglequadGauss7w = @SVector [ 0.225000000000000, 0.125939180544827, 0.125939180544827, @@ -145,7 +149,7 @@ const Trianglegauss7w = [ 0.132394152788506, ] -const Trianglegauss12a = [ +const trianglequadGauss12a = @SVector [ 0.873821971016996, 0.063089014491502, 0.063089014491502, @@ -160,7 +164,7 @@ const Trianglegauss12a = [ 0.053145049844816, ] -const Trianglegauss12b = [ +const trianglequadGauss12b = @SVector [ 0.063089014491502, 0.873821971016996, 0.063089014491502, @@ -175,7 +179,7 @@ const Trianglegauss12b = [ 0.310352451033785, ] -const Trianglegauss12c = [ +const trianglequadGauss12c = @SVector [ 0.063089014491502, 0.063089014491502, 0.873821971016996, @@ -190,7 +194,7 @@ const Trianglegauss12c = [ 0.636502499121399, ] -const Trianglegauss12w = [ +const trianglequadGauss12w = @SVector [ 0.050844906370207, 0.050844906370207, 0.050844906370207, @@ -205,7 +209,7 @@ const Trianglegauss12w = [ 0.082851075618374, ] -const Trianglegauss13a = [ +const trianglequadGauss13a = @SVector [ 0.333333333333333, 0.479308067841923, 0.260345966079038, @@ -221,7 +225,7 @@ const Trianglegauss13a = [ 0.0486903154253160, ] -const Trianglegauss13b = [ +const trianglequadGauss13b = @SVector [ 0.333333333333333, 0.260345966079038, 0.479308067841923, @@ -237,7 +241,7 @@ const Trianglegauss13b = [ 0.312865496004875, ] -const Trianglegauss13c = [ +const trianglequadGauss13c = @SVector [ 0.333333333333333, 0.260345966079038, 0.260345966079038, @@ -253,7 +257,7 @@ const Trianglegauss13c = [ 0.638444188569809, ] -const Trianglegauss13w = [ +const trianglequadGauss13w = @SVector [ -0.149570044467670, 0.175615257433204, 0.175615257433204, @@ -269,7 +273,7 @@ const Trianglegauss13w = [ 0.077113760890257, ] -const Trianglegauss36a = [ +const trianglequadGauss36a = @SVector [ 0.0242935351590, 0.0265193427722, 0.9492126023551, @@ -308,7 +312,7 @@ const Trianglegauss36a = [ 0.2305424298836 ] -const Trianglegauss36b = [ +const trianglequadGauss36b = @SVector [ 0.9493059293846, 0.0242695130640, 0.0265067966437, @@ -347,7 +351,7 @@ const Trianglegauss36b = [ 0.3456013949376 ] -const Trianglegauss36w = [ +const trianglequadGauss36w = (@SVector [ 0.0166240998757, 0.0166811699778, 0.0166830569067, @@ -384,9 +388,9 @@ const Trianglegauss36w = [ 0.1025573374896, 0.1033159661413, 0.1035854367193 -] / 2 +]) ./ 2 -const Trianglegauss78a = [ +const trianglequadGauss78a = @SVector [ 0.0089411337112, 0.9792622629807, 0.0105475382112, @@ -467,7 +471,7 @@ const Trianglegauss78a = [ 0.3729077987144 ] -const Trianglegauss78b = [ +const trianglequadGauss78b = @SVector [ 0.0086983293702, 0.0102644133744, 0.9785514202515, @@ -548,7 +552,7 @@ const Trianglegauss78b = [ 0.3753750277549 ] -const Trianglegauss78w = [ +const trianglequadGauss78w = (@SVector [ 0.0021744545399, 0.0028987135265, 0.0030846029337, @@ -627,9 +631,9 @@ const Trianglegauss78w = [ 0.0634133183449, 0.0635311861108, 0.0637206605672 -] / 2 +]) ./ 2 -const Trianglegauss105a = [ +const trianglequadGauss105a = @SVector [ 0.0087809303836, 0.9903675314220, 0.0027029276450, @@ -737,7 +741,7 @@ const Trianglegauss105a = [ 0.4098894602340 ] -const Trianglegauss105b = [ +const trianglequadGauss105b = @SVector [ 0.9903676436772, 0.0087809216232, 0.0335914404439, @@ -845,7 +849,7 @@ const Trianglegauss105b = [ 0.4098894317792, ] -const Trianglegauss105w = [ +const trianglequadGauss105w = (@SVector [ 0.0006438298261, 0.0006438413076, 0.0010134735710, @@ -951,9 +955,9 @@ const Trianglegauss105w = [ 0.0392705643548, 0.0392705802517, 0.0398766879831 -] / 2 +]) ./ 2 -const Trianglegauss120a = [ +const trianglequadGauss120a = @SVector [ 0.0082881595033, 0.4618422030241, 0.0071066441239, @@ -1076,7 +1080,7 @@ const Trianglegauss120a = [ 0.3361523347440 ] -const Trianglegauss120b = [ +const trianglequadGauss120b = @SVector [ 0.9848202768869, 0.5381577969759, 0.0080842361390, @@ -1199,7 +1203,7 @@ const Trianglegauss120b = [ 0.2778500044356 ] -const Trianglegauss120w = [ +const trianglequadGauss120w = (@SVector [ 0.0014873417859, 0.0014889035262, 0.0015005944380, @@ -1320,9 +1324,9 @@ const Trianglegauss120w = [ 0.0364656225016, 0.0365172708706, 0.0371924811018 -] / 2 +]) ./ 2 -const Trianglegauss400a = [ +const trianglequadGauss400a = @SVector [ 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, @@ -1404,7 +1408,7 @@ const Trianglegauss400a = [ 9.9656429959254744e-01, 9.9656429959254744e-01, 9.9656429959254744e-01, 9.9656429959254744e-01, 9.9656429959254744e-01, 9.9656429959254744e-01, 9.9656429959254744e-01, 9.9656429959254744e-01, 9.9656429959254744e-01, 9.9656429959254744e-01 ] -const Trianglegauss400b = [ +const trianglequadGauss400b = @SVector [ 3.4238963701627881e-03, 1.7952145528977594e-02, 4.3732017769028415e-02, 8.0165141146159247e-02, 1.2639828298375833e-01, 1.8134795437803311e-01, 2.4372624180423841e-01, 3.1207107651919747e-01, 3.8478052619624081e-01, 4.6015035032944468e-01, 5.3641394926310271e-01, 6.1178377339630663e-01, 6.8449322307334992e-01, 7.5283805778830903e-01, 8.1521634521451436e-01, @@ -1486,7 +1490,7 @@ const Trianglegauss400b = [ 1.8493113036459207e-03, 2.1091521745159178e-03, 2.3598212843598133e-03, 2.5954431871045409e-03, 2.8104951487431724e-03, 2.9999366212862710e-03, 3.1593274647212558e-03, 3.2849323021439448e-03, 3.3738095753870568e-03, 3.4238963701627881e-03 ] -const Trianglegauss400w = [ +const trianglequadGauss400w = (@SVector [ 3.0918731028920757e-04, 7.1269681990617786e-04, 1.1001132168084852e-03, 1.4617975077641802e-03, 1.7892294090025515e-03, 2.0747266161171671e-03, 2.3115952886863795e-03, 2.4942827316709218e-03, 2.6185066288951639e-03, 2.6813551585055523e-03, 2.6813551585055523e-03, 2.6185066288951639e-03, 2.4942827316709218e-03, 2.3115952886863795e-03, 2.0747266161171671e-03, @@ -1566,9 +1570,9 @@ const Trianglegauss400w = [ 1.0659372088425298e-06, 2.4570544575427101e-06, 3.7926899737209111e-06, 5.0396128931087845e-06, 6.1684491528038026e-06, 7.1527136615879260e-06, 7.9693291024496169e-06, 8.5991523085940462e-06, 9.0274198017032943e-06, 9.2440929444984856e-06, 9.2440929444984856e-06, 9.0274198017032943e-06, 8.5991523085940462e-06, 7.9693291024496169e-06, 7.1527136615879260e-06, - 6.1684491528038026e-06, 5.0396128931087845e-06, 3.7926899737209111e-06, 2.4570544575427101e-06, 1.0659372088425298e-06 ] / 2 + 6.1684491528038026e-06, 5.0396128931087845e-06, 3.7926899737209111e-06, 2.4570544575427101e-06, 1.0659372088425298e-06 ]) ./ 2 -const Trianglegauss900a = [ +const trianglequadGauss900a = @SVector [ 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, @@ -1750,7 +1754,7 @@ const Trianglegauss900a = [ 9.9844674203732475e-01, 9.9844674203732475e-01, 9.9844674203732475e-01, 9.9844674203732475e-01, 9.9844674203732475e-01, 9.9844674203732475e-01, 9.9844674203732475e-01, 9.9844674203732475e-01, 9.9844674203732475e-01, 9.9844674203732475e-01 ] -const Trianglegauss900b = [ +const trianglequadGauss900b = @SVector [ 1.5508453523766334e-03, 8.1532545513458305e-03, 1.9958019337560785e-02, 3.6842661103375078e-02, 5.8628525212516977e-02, 8.5084754640270099e-02, 1.1593093317123887e-01, 1.5084009454090175e-01, 1.8944219814142918e-01, 2.3132805481318516e-01, 2.7605366563067762e-01, 3.2314492896492530e-01, 3.7210266635129818e-01, 4.2240791404863498e-01, 4.7352742426565486e-01, @@ -1932,7 +1936,7 @@ const Trianglegauss900b = [ 1.1933868468703306e-03, 1.2585475987258832e-03, 1.3185999002349133e-03, 1.3729071863272115e-03, 1.4208937943295222e-03, 1.4620510712185931e-03, 1.4959427806874798e-03, 1.5222097843898060e-03, 1.5405741538946657e-03, 1.5508453523766334e-03 ] -const Trianglegauss900w = [ +const trianglequadGauss900w = (@SVector [ 6.3393472058975955e-05, 1.4691582105612550e-04, 2.2900583486229500e-04, 3.0867923073618977e-04, 3.8508275109363019e-04, 4.5740496279664960e-04, 5.2487882297804151e-04, 5.8678895194472442e-04, 6.4247903025549066e-04, 6.9135870330888595e-04, 7.3290982143010656e-04, 7.6669192547512504e-04, 7.9234691283714938e-04, 8.0960283204461417e-04, 8.1827676484210919e-04, @@ -2112,66 +2116,52 @@ const Trianglegauss900w = [ 1.1401689926257982e-06, 1.1927229445739803e-06, 1.2326337497523175e-06, 1.2594783402385204e-06, 1.2729721547988522e-06, 1.2729721547988522e-06, 1.2594783402385204e-06, 1.2326337497523175e-06, 1.1927229445739803e-06, 1.1401689926257982e-06, 1.0755289849393060e-06, 9.9948813249661188e-07, 9.1285250744394168e-07, 8.1654050927813780e-07, 7.1157315730367923e-07, - 5.9906334934258464e-07, 4.8020435429046972e-07, 3.5625849784750484e-07, 2.2855317093104202e-07, 9.8619596931444806e-08 ] / 2 - -const triangleGaussA = Vector{Float64}[ - Trianglegauss1a, - Trianglegauss3a, - Trianglegauss4a, - Trianglegauss6a, - Trianglegauss7a, - Trianglegauss12a, - Trianglegauss13a, - Trianglegauss36a, - Trianglegauss78a, - Trianglegauss105a, - Trianglegauss120a, - Trianglegauss400a, - Trianglegauss900a -] + 5.9906334934258464e-07, 4.8020435429046972e-07, 3.5625849784750484e-07, 2.2855317093104202e-07, 9.8619596931444806e-08 ]) ./ 2 -const triangleGaussB = Vector{Float64}[ - Trianglegauss1b, - Trianglegauss3b, - Trianglegauss4b, - Trianglegauss6b, - Trianglegauss7b, - Trianglegauss12b, - Trianglegauss13b, - Trianglegauss36b, - Trianglegauss78b, - Trianglegauss105b, - Trianglegauss120b, - Trianglegauss400b, - Trianglegauss900b +const trianglequadGaussA = @SVector [ + trianglequadGauss1a, + trianglequadGauss3a, + trianglequadGauss4a, + trianglequadGauss6a, + trianglequadGauss7a, + trianglequadGauss12a, + trianglequadGauss13a, + trianglequadGauss36a, + trianglequadGauss78a, + trianglequadGauss105a, + trianglequadGauss120a, + trianglequadGauss400a, + trianglequadGauss900a ] -const triangleGaussW = Vector{Float64}[ - Trianglegauss1w, - Trianglegauss3w, - Trianglegauss4w, - Trianglegauss6w, - Trianglegauss7w, - Trianglegauss12w, - Trianglegauss13w, - Trianglegauss36w, - Trianglegauss78w, - Trianglegauss105w, - Trianglegauss120w, - Trianglegauss400w, - Trianglegauss900w +const trianglequadGaussB = @SVector [ + trianglequadGauss1b, + trianglequadGauss3b, + trianglequadGauss4b, + trianglequadGauss6b, + trianglequadGauss7b, + trianglequadGauss12b, + trianglequadGauss13b, + trianglequadGauss36b, + trianglequadGauss78b, + trianglequadGauss105b, + trianglequadGauss120b, + trianglequadGauss400b, + trianglequadGauss900b ] -""" - trgauss(n) -> (u,w) - -Returns the n-th triangle quadrature rule. Returns a Matrix u of size (Q,2) -with Q the number of quadrature points and a Vector w of size (Q,) containing -the quadrature weights. -""" -function trgauss(n) - @assert 1 <= n <= length(triangleGaussA) - u = copy(transpose([triangleGaussA[n] triangleGaussB[n]])) - w = triangleGaussW[n] - return u, w/2 -end +const trianglequadGaussW = @SVector [ + trianglequadGauss1w, + trianglequadGauss3w, + trianglequadGauss4w, + trianglequadGauss6w, + trianglequadGauss7w, + trianglequadGauss12w, + trianglequadGauss13w, + trianglequadGauss36w, + trianglequadGauss78w, + trianglequadGauss105w, + trianglequadGauss120w, + trianglequadGauss400w, + trianglequadGauss900w +] \ No newline at end of file diff --git a/src/quadrature/triangle/TriangleQuadratures.jl b/src/quadrature/triangle/TriangleQuadratures.jl new file mode 100644 index 0000000..0d338e3 --- /dev/null +++ b/src/quadrature/triangle/TriangleQuadratures.jl @@ -0,0 +1,43 @@ +include("TriangleDunavant.jl") +include("TriangleGauss.jl") + +export trgauss +export TriangleQuadDunavant +export TriangleQuadGauss + +abstract type TriangleQuadrature end + +struct TriangleQuadDunavant{I} <: TriangleQuadrature + order::I +end + +struct TriangleQuadGauss{I} <: TriangleQuadrature + order::I +end + +""" + trgauss(n) -> (u,w) + +Returns the n-th triangle quadrature rule. Returns a Matrix u of size (Q,2) +with Q the number of quadrature points and a Vector w of size (Q,) containing +the quadrature weights. +""" +function trgauss(n) + return trgauss(TriangleQuadGauss(n)) +end + +function trgauss(trqd::TriangleQuadDunavant) + n = trqd.order + @assert 1 <= n <= length(trianglequadDunavantA) + u = copy(transpose([trianglequadDunavantA[n] trianglequadDunavantB[n]])) + w = trianglequadDunavantW[n] + return u, w/2 +end + +function trgauss(trqd::TriangleQuadGauss) + n = trqd.order + @assert 1 <= n <= length(trianglequadGaussA) + u = copy(transpose([trianglequadGaussA[n] trianglequadGaussB[n]])) + w = trianglequadGaussW[n] + return u, w/2 +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 5dc5917..cce028d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -39,6 +39,7 @@ include("test_isinclosure.jl") include("test_permute_vertices.jl") include("test_trgauss.jl") +include("test_trdunavant.jl") include("test_chartquad.jl") include("test_spherequad.jl") include("test_nbd.jl") diff --git a/test/test_trdunavant.jl b/test/test_trdunavant.jl new file mode 100644 index 0000000..49962eb --- /dev/null +++ b/test/test_trdunavant.jl @@ -0,0 +1,91 @@ +using CompScienceMeshes +using Test +using StaticArrays + +#= The integrand used to test the quadrature is given by + +f(x,y) = x^o + y^o + x^(floor(o/2)) * y^(ceil(o/2)), + +where o is the order of the polynomial is used for testing. + +Each quadrature rule 'r' is supposed to integrate +the corresponding polynomial resulting by 'o:=r', exactly. + +Thus, the functions used are + +functions = [ + "x¹ + y¹ + y", + "x² + y² + xy", + "x³ + y³ + xy²", + "x⁴ + y⁴ + x²y²", + "x⁵ + y⁵ + x²y³", + "x⁶ + y⁶ + x³y³", + "x⁷ + y⁷ + x³y⁴", + "x⁸ + y⁸ + x⁴y⁴", + "x⁹ + y⁹ + x⁴y⁵", + "x¹⁰ + y¹⁰ + x⁵y⁵", + "x¹¹ + y¹¹ + x⁵y⁶", + "x¹² + y¹² + x⁶y⁶", + "x¹³ + y¹³ + x⁶y⁷", + "x¹⁴ + y¹⁴ + x⁷y⁷", + "x¹⁵ + y¹⁵ + x⁷y⁸", + "x¹⁶ + y¹⁶ + x⁸y⁸", + "x¹⁷ + y¹⁷ + x⁸y⁹", + "x¹⁸ + y¹⁸ + x⁹y⁹", + "x¹⁹ + y¹⁹ + x⁹y¹⁰", + "x²⁰ + y²⁰ + x¹⁰y¹⁰", +] +=# + +ignd(x, y, order::Int) = (x^order + y^order + x^(ceil(order/2)) * y^(floor(order/2))) + +N = length(CompScienceMeshes.trianglequadDunavantW) + +for T in @SVector [Float32, Float64] + + local p1 = T.([0.0, 0.0]) + local p2 = T.([1.0, 0.0]) + local p3 = T.([0.0, 1.0]) + + # Analytical Results + local ana = @SVector (T[ + 0.5, # 1 / 2, + 0.2083333333333333400, # 5 / 24, + 0.1166666666666666700, # 7 / 60, + 0.0722222222222222200, # 13 / 180, + 0.05, # 1 / 20, + 0.0366071428571428600, # 41 / 1120, + 0.0281746031746031750, # 71 / 2520, + 0.0223809523809523800, # 47 / 2100, + 0.0182539682539682550, # 23 / 1260, + 0.0151815776815776810, # 505 / 33264, + 0.0128343878343878340, # 925 / 72072, + 0.0109949574235288520, # 1849 / 168168, + 0.0095265845265845270, # 3433 / 360360, + 0.0083345473970473980, # 1373 / 164736, + 0.0073535125005713240, # 12871 / 1750320, + 0.0065362016342408500, # 25741 / 3938220, + 0.0058480734951323185, # 8314191 / 1421697420, + 0.0052632120201779640, # 97241 / 18475600, + 0.0047619305359243440, # 3879897 / 814773960, + 0.0043290160444677760, # 123171 / 28452424, + ]) + + J = zeros(N) + for _n in 1 : N + u, w = trgauss(TriangleQuadDunavant(_n)) + for _g in eachindex(w) + _p = u[1,_g]*p1 + u[2,_g]*p2 + (1-u[1,_g]-u[2,_g])*p3 + x, y = _p[1], _p[2] + J[_n] += w[_g] * ignd(x, y, _n) + end + end + + # Tolerance factor of 200 for Float64 necessary + T == Float32 ? tol = T(1) : tol = T(200) + + for _q in 1 : N + #@show abs(J[_q]-ana[_q])/abs(ana[_q]) + @test abs(J[_q]-ana[_q])/abs(ana[_q]) <= tol*eps(T) + end +end \ No newline at end of file diff --git a/test/test_trgauss.jl b/test/test_trgauss.jl index 8f3cd7e..f927fc7 100644 --- a/test/test_trgauss.jl +++ b/test/test_trgauss.jl @@ -9,7 +9,7 @@ for T in [Float32, Float64] local p2 = T.([1.0, 0.0]) local p3 = T.([0.0, 1.0]) - N = length(CompScienceMeshes.triangleGaussW) + N = length(CompScienceMeshes.trianglequadGaussW) J = zeros(N) for _n in 1 : N u, w = trgauss(_n) From 8b99012a632715f41ddfbebf8a1eea837ef54a92 Mon Sep 17 00:00:00 2001 From: Kristof Cools Date: Wed, 28 Aug 2024 11:36:14 +0200 Subject: [PATCH 2/3] Name change Gaussian to Legacy and use testitems --- src/quadrature/triangle/TriangleQuadratures.jl | 8 ++++---- test/test_trdunavant.jl | 3 +++ test/test_trgauss.jl | 6 ++++-- 3 files changed, 11 insertions(+), 6 deletions(-) diff --git a/src/quadrature/triangle/TriangleQuadratures.jl b/src/quadrature/triangle/TriangleQuadratures.jl index 0d338e3..b6f564d 100644 --- a/src/quadrature/triangle/TriangleQuadratures.jl +++ b/src/quadrature/triangle/TriangleQuadratures.jl @@ -3,7 +3,7 @@ include("TriangleGauss.jl") export trgauss export TriangleQuadDunavant -export TriangleQuadGauss +export TriangleQuadLegacy abstract type TriangleQuadrature end @@ -11,7 +11,7 @@ struct TriangleQuadDunavant{I} <: TriangleQuadrature order::I end -struct TriangleQuadGauss{I} <: TriangleQuadrature +struct TriangleQuadLegacy{I} <: TriangleQuadrature order::I end @@ -23,7 +23,7 @@ with Q the number of quadrature points and a Vector w of size (Q,) containing the quadrature weights. """ function trgauss(n) - return trgauss(TriangleQuadGauss(n)) + return trgauss(TriangleQuadLegacy(n)) end function trgauss(trqd::TriangleQuadDunavant) @@ -34,7 +34,7 @@ function trgauss(trqd::TriangleQuadDunavant) return u, w/2 end -function trgauss(trqd::TriangleQuadGauss) +function trgauss(trqd::TriangleQuadLegacy) n = trqd.order @assert 1 <= n <= length(trianglequadGaussA) u = copy(transpose([trianglequadGaussA[n] trianglequadGaussB[n]])) diff --git a/test/test_trdunavant.jl b/test/test_trdunavant.jl index 49962eb..ed97967 100644 --- a/test/test_trdunavant.jl +++ b/test/test_trdunavant.jl @@ -37,6 +37,8 @@ functions = [ ] =# +@testitem "accuracy" begin +using StaticArrays ignd(x, y, order::Int) = (x^order + y^order + x^(ceil(order/2)) * y^(floor(order/2))) N = length(CompScienceMeshes.trianglequadDunavantW) @@ -88,4 +90,5 @@ for T in @SVector [Float32, Float64] #@show abs(J[_q]-ana[_q])/abs(ana[_q]) @test abs(J[_q]-ana[_q])/abs(ana[_q]) <= tol*eps(T) end +end end \ No newline at end of file diff --git a/test/test_trgauss.jl b/test/test_trgauss.jl index f927fc7..b324c7a 100644 --- a/test/test_trgauss.jl +++ b/test/test_trgauss.jl @@ -1,6 +1,7 @@ -using CompScienceMeshes -using Test +# using CompScienceMeshes +# using Test +@testitem "accuracy" begin ignd(x,y) = x^2 * cos(x*y) for T in [Float32, Float64] Iref = 0.082739596439277 @@ -25,4 +26,5 @@ for T in [Float32, Float64] for _q in 1 : N @test abs(J[_q]-Iref)/abs(Iref) <= E[_q]*2.0 end +end end \ No newline at end of file From f388b1e0bfda05b810ddf3f22242394f8bf47f79 Mon Sep 17 00:00:00 2001 From: Kristof Cools Date: Wed, 28 Aug 2024 14:02:04 +0200 Subject: [PATCH 3/3] Revert to using Array to store quadrules --- src/quadrature/triangle/TriangleDunavant.jl | 126 +++++++++--------- src/quadrature/triangle/TriangleGauss.jl | 98 +++++++------- .../triangle/TriangleQuadratures.jl | 18 ++- 3 files changed, 128 insertions(+), 114 deletions(-) diff --git a/src/quadrature/triangle/TriangleDunavant.jl b/src/quadrature/triangle/TriangleDunavant.jl index 71244d6..e827559 100644 --- a/src/quadrature/triangle/TriangleDunavant.jl +++ b/src/quadrature/triangle/TriangleDunavant.jl @@ -3,58 +3,58 @@ # Rules for the Triangle.” International Journal for Numerical Methods # in Engineering 21, no. 6 (1985): 1129–48. https://doi.org/10.1002/nme.1620210612. -const trianglequadDunavant1a = @SVector [ +const trianglequadDunavant1a = [ 0.333333333333333, ] -const trianglequadDunavant1b = @SVector [ +const trianglequadDunavant1b = [ 0.333333333333333, ] -const trianglequadDunavant1w = @SVector [ +const trianglequadDunavant1w = [ 1.000000000000000, ] -const trianglequadDunavant3a = @SVector [ +const trianglequadDunavant3a = [ 0.666666666666667, 0.166666666666667, 0.166666666666667, ] -const trianglequadDunavant3b = @SVector [ +const trianglequadDunavant3b = [ 0.166666666666667, 0.166666666666667, 0.666666666666667, ] -const trianglequadDunavant3w = @SVector [ +const trianglequadDunavant3w = [ 0.333333333333333, 0.333333333333333, 0.333333333333333, ] -const trianglequadDunavant4a = @SVector [ +const trianglequadDunavant4a = [ 0.333333333333333, 0.600000000000000, 0.200000000000000, 0.200000000000000, ] -const trianglequadDunavant4b = @SVector [ +const trianglequadDunavant4b = [ 0.333333333333333, 0.200000000000000, 0.200000000000000, 0.600000000000000, ] -const trianglequadDunavant4w = @SVector [ +const trianglequadDunavant4w = [ -0.562500000000000, 0.520833333333333, 0.520833333333333, 0.520833333333333, ] -const trianglequadDunavant6a = @SVector [ +const trianglequadDunavant6a = [ 0.108103018168070, 0.445948490915965, 0.445948490915965, @@ -63,7 +63,7 @@ const trianglequadDunavant6a = @SVector [ 0.091576213509771, ] -const trianglequadDunavant6b = @SVector [ +const trianglequadDunavant6b = [ 0.445948490915965, 0.445948490915965, 0.108103018168070, @@ -72,7 +72,7 @@ const trianglequadDunavant6b = @SVector [ 0.816847572980459, ] -const trianglequadDunavant6w = @SVector [ +const trianglequadDunavant6w = [ 0.223381589678011, 0.223381589678011, 0.223381589678011, @@ -81,7 +81,7 @@ const trianglequadDunavant6w = @SVector [ 0.109951743655322, ] -const trianglequadDunavant7a = @SVector [ +const trianglequadDunavant7a = [ 0.333333333333333, 0.059715871789770, 0.470142064105115, @@ -91,7 +91,7 @@ const trianglequadDunavant7a = @SVector [ 0.101286507323456, ] -const trianglequadDunavant7b = @SVector [ +const trianglequadDunavant7b = [ 0.333333333333333, 0.470142064105115, 0.470142064105115, @@ -101,7 +101,7 @@ const trianglequadDunavant7b = @SVector [ 0.797426985353087, ] -const trianglequadDunavant7w = @SVector [ +const trianglequadDunavant7w = [ 0.225000000000000, 0.132394152788506, 0.132394152788506, @@ -111,7 +111,7 @@ const trianglequadDunavant7w = @SVector [ 0.125939180544827, ] -const trianglequadDunavant12a = @SVector [ +const trianglequadDunavant12a = [ 0.501426509658179, 0.249286745170910, 0.249286745170910, @@ -126,7 +126,7 @@ const trianglequadDunavant12a = @SVector [ 0.053145049844817, ] -const trianglequadDunavant12b = @SVector [ +const trianglequadDunavant12b = [ 0.249286745170910, 0.249286745170910, 0.501426509658179, @@ -141,7 +141,7 @@ const trianglequadDunavant12b = @SVector [ 0.636502499121399, ] -const trianglequadDunavant12w = @SVector [ +const trianglequadDunavant12w = [ 0.116786275726379, 0.116786275726379, 0.116786275726379, @@ -156,7 +156,7 @@ const trianglequadDunavant12w = @SVector [ 0.082851075618374, ] -const trianglequadDunavant13a = @SVector [ +const trianglequadDunavant13a = [ 0.333333333333333, 0.479308067841920, 0.260345966079040, @@ -172,7 +172,7 @@ const trianglequadDunavant13a = @SVector [ 0.048690315425316, ] -const trianglequadDunavant13b = @SVector [ +const trianglequadDunavant13b = [ 0.333333333333333, 0.260345966079040, 0.260345966079040, @@ -188,7 +188,7 @@ const trianglequadDunavant13b = @SVector [ 0.638444188569810, ] -const trianglequadDunavant13w = @SVector [ +const trianglequadDunavant13w = [ -0.149570044467682, 0.175615257433208, 0.175615257433208, @@ -204,7 +204,7 @@ const trianglequadDunavant13w = @SVector [ 0.077113760890257, ] -const trianglequadDunavant16a = @SVector [ +const trianglequadDunavant16a = [ 0.333333333333333, 0.081414823414554, 0.459292588292723, @@ -223,7 +223,7 @@ const trianglequadDunavant16a = @SVector [ 0.008394777409958, ] -const trianglequadDunavant16b = @SVector [ +const trianglequadDunavant16b = [ 0.333333333333333, 0.459292588292723, 0.459292588292723, @@ -242,7 +242,7 @@ const trianglequadDunavant16b = @SVector [ 0.728492392955404, ] -const trianglequadDunavant16w = @SVector [ +const trianglequadDunavant16w = [ 0.144315607677787, 0.095091634267285, 0.095091634267285, @@ -261,7 +261,7 @@ const trianglequadDunavant16w = @SVector [ 0.027230314174435, ] -const trianglequadDunavant19a = @SVector [ +const trianglequadDunavant19a = [ 0.333333333333333, 0.020634961602525, 0.489682519198738, @@ -283,7 +283,7 @@ const trianglequadDunavant19a = @SVector [ 0.036838412054736, ] -const trianglequadDunavant19b = @SVector [ +const trianglequadDunavant19b = [ 0.333333333333333, 0.489682519198738, 0.489682519198738, @@ -305,7 +305,7 @@ const trianglequadDunavant19b = @SVector [ 0.741198598784498, ] -const trianglequadDunavant19w = @SVector [ +const trianglequadDunavant19w = [ 0.097135796282799, 0.031334700227139, 0.031334700227139, @@ -327,7 +327,7 @@ const trianglequadDunavant19w = @SVector [ 0.043283539377289, ] -const trianglequadDunavant25a = @SVector [ +const trianglequadDunavant25a = [ 0.333333333333333, 0.028844733232685, 0.485577633383657, @@ -355,7 +355,7 @@ const trianglequadDunavant25a = @SVector [ 0.009540815400299, ] -const trianglequadDunavant25b = @SVector [ +const trianglequadDunavant25b = [ 0.333333333333333, 0.485577633383657, 0.485577633383657, @@ -383,7 +383,7 @@ const trianglequadDunavant25b = @SVector [ 0.923655933587500, ] -const trianglequadDunavant25w = @SVector [ +const trianglequadDunavant25w = [ 0.090817990382754, 0.036725957756467, 0.036725957756467, @@ -411,7 +411,7 @@ const trianglequadDunavant25w = @SVector [ 0.009421666963733, ] -const trianglequadDunavant27a = @SVector [ +const trianglequadDunavant27a = [ -0.069222096541517, 0.534611048270758, 0.534611048270758, @@ -441,7 +441,7 @@ const trianglequadDunavant27a = @SVector [ 0.021022016536166, ] -const trianglequadDunavant27b = @SVector [ +const trianglequadDunavant27b = [ 0.534611048270758, 0.534611048270758, -0.069222096541517, @@ -471,7 +471,7 @@ const trianglequadDunavant27b = @SVector [ 0.807489003159792, ] -const trianglequadDunavant27w = @SVector [ +const trianglequadDunavant27w = [ 0.000927006328961, 0.000927006328961, 0.000927006328961, @@ -501,7 +501,7 @@ const trianglequadDunavant27w = @SVector [ 0.020707659639141, ] -const trianglequadDunavant33a = @SVector [ +const trianglequadDunavant33a = [ 0.023565220452390, 0.488217389773805, 0.488217389773805, @@ -537,7 +537,7 @@ const trianglequadDunavant33a = @SVector [ 0.025734050548330, ] -const trianglequadDunavant33b = @SVector [ +const trianglequadDunavant33b = [ 0.488217389773805, 0.488217389773805, 0.023565220452390, @@ -573,7 +573,7 @@ const trianglequadDunavant33b = @SVector [ 0.858014033544073, ] -const trianglequadDunavant33w = @SVector [ +const trianglequadDunavant33w = [ 0.025731066440455, 0.025731066440455, 0.025731066440455, @@ -609,7 +609,7 @@ const trianglequadDunavant33w = @SVector [ 0.017316231108659, ] -const trianglequadDunavant37a = @SVector [ +const trianglequadDunavant37a = [ 0.333333333333333, 0.009903630120591, 0.495048184939705, @@ -649,7 +649,7 @@ const trianglequadDunavant37a = @SVector [ 0.022233076674090, ] -const trianglequadDunavant37b = @SVector [ +const trianglequadDunavant37b = [ 0.333333333333333, 0.495048184939705, 0.495048184939705, @@ -689,7 +689,7 @@ const trianglequadDunavant37b = @SVector [ 0.851409537834241, ] -const trianglequadDunavant37w = @SVector [ +const trianglequadDunavant37w = [ 0.052520923400802, 0.011280145209330, 0.011280145209330, @@ -729,7 +729,7 @@ const trianglequadDunavant37w = @SVector [ 0.015521786839045, ] -const trianglequadDunavant42a = @SVector [ +const trianglequadDunavant42a = [ 0.022072179275643, 0.488963910362179, 0.488963910362179, @@ -774,7 +774,7 @@ const trianglequadDunavant42a = @SVector [ 0.001268330932872, ] -const trianglequadDunavant42b = @SVector [ +const trianglequadDunavant42b = [ 0.488963910362179, 0.488963910362179, 0.022072179275643, @@ -819,7 +819,7 @@ const trianglequadDunavant42b = @SVector [ 0.879757171370171, ] -const trianglequadDunavant42w = @SVector [ +const trianglequadDunavant42w = [ 0.021883581369429, 0.021883581369429, 0.021883581369429, @@ -864,7 +864,7 @@ const trianglequadDunavant42w = @SVector [ 0.005010228838501, ] -const trianglequadDunavant48a = @SVector [ +const trianglequadDunavant48a = [ -0.013945833716486, 0.506972916858243, 0.506972916858243, @@ -915,7 +915,7 @@ const trianglequadDunavant48a = @SVector [ 0.012459809331199, ] -const trianglequadDunavant48b = @SVector [ +const trianglequadDunavant48b = [ 0.506972916858243, 0.506972916858243, -0.013945833716486, @@ -966,7 +966,7 @@ const trianglequadDunavant48b = @SVector [ 0.883964574092416, ] -const trianglequadDunavant48w = @SVector [ +const trianglequadDunavant48w = [ 0.001916875642849, 0.001916875642849, 0.001916875642849, @@ -1017,7 +1017,7 @@ const trianglequadDunavant48w = @SVector [ 0.007673942631049, ] -const trianglequadDunavant52a = @SVector [ +const trianglequadDunavant52a = [ 0.333333333333333, 0.005238916103123, 0.497380541948438, @@ -1072,7 +1072,7 @@ const trianglequadDunavant52a = @SVector [ 0.014317320230681, ] -const trianglequadDunavant52b = @SVector [ +const trianglequadDunavant52b = [ 0.333333333333333, 0.497380541948438, 0.497380541948438, @@ -1127,7 +1127,7 @@ const trianglequadDunavant52b = @SVector [ 0.900399064086661, ] -const trianglequadDunavant52w = @SVector [ +const trianglequadDunavant52w = [ 0.046875697427642, 0.006405878578585, 0.006405878578585, @@ -1182,7 +1182,7 @@ const trianglequadDunavant52w = @SVector [ 0.006850054546542, ] -const trianglequadDunavant61a = @SVector [ +const trianglequadDunavant61a = [ 0.333333333333333, 0.005658918886452, 0.497170540556774, @@ -1246,7 +1246,7 @@ const trianglequadDunavant61a = @SVector [ 0.014663182224828, ] -const trianglequadDunavant61b = @SVector [ +const trianglequadDunavant61b = [ 0.333333333333333, 0.497170540556774, 0.497170540556774, @@ -1310,7 +1310,7 @@ const trianglequadDunavant61b = @SVector [ 0.904625504095608, ] -const trianglequadDunavant61w = @SVector [ +const trianglequadDunavant61w = [ 0.033437199290803, 0.005093415440507, 0.005093415440507, @@ -1374,7 +1374,7 @@ const trianglequadDunavant61w = @SVector [ 0.006665632004165, ] -const trianglequadDunavant70a = @SVector [ +const trianglequadDunavant70a = [ 0.333333333333333, 0.013310382738157, 0.493344808630921, @@ -1447,7 +1447,7 @@ const trianglequadDunavant70a = @SVector [ -0.035222015287949, ] -const trianglequadDunavant70b = @SVector [ +const trianglequadDunavant70b = [ 0.333333333333333, 0.493344808630921, 0.493344808630921, @@ -1520,7 +1520,7 @@ const trianglequadDunavant70b = @SVector [ 1.014347260005363, ] -const trianglequadDunavant70w = @SVector [ +const trianglequadDunavant70w = [ 0.030809939937647, 0.009072436679404, 0.009072436679404, @@ -1593,7 +1593,7 @@ const trianglequadDunavant70w = @SVector [ 0.000046187660794, ] -const trianglequadDunavant73a = @SVector [ +const trianglequadDunavant73a = [ 0.333333333333333, 0.020780025853987, 0.489609987073006, @@ -1669,7 +1669,7 @@ const trianglequadDunavant73a = @SVector [ 0.010161119296278, ] -const trianglequadDunavant73b = @SVector [ +const trianglequadDunavant73b = [ 0.333333333333333, 0.489609987073006, 0.489609987073006, @@ -1745,7 +1745,7 @@ const trianglequadDunavant73b = @SVector [ 0.924344252620784, ] -const trianglequadDunavant73w = @SVector [ +const trianglequadDunavant73w = [ 0.032906331388919, 0.010330731891272, 0.010330731891272, @@ -1821,7 +1821,7 @@ const trianglequadDunavant73w = @SVector [ 0.003799928855302, ] -const trianglequadDunavant79a = @SVector [ +const trianglequadDunavant79a = [ 0.333333333333333, -0.001900928704400, 0.500950464352200, @@ -1903,7 +1903,7 @@ const trianglequadDunavant79a = @SVector [ 0.010547719294141, ] -const trianglequadDunavant79b = @SVector [ +const trianglequadDunavant79b = [ 0.333333333333333, 0.500950464352200, 0.500950464352200, @@ -1985,7 +1985,7 @@ const trianglequadDunavant79b = @SVector [ 0.929756171556853, ] -const trianglequadDunavant79w = @SVector [ +const trianglequadDunavant79w = [ 0.033057055541624, 0.000867019185663, 0.000867019185663, @@ -2067,7 +2067,7 @@ const trianglequadDunavant79w = @SVector [ 0.003573909385950, ] -const trianglequadDunavantA = @SVector [ +const trianglequadDunavantA = [ trianglequadDunavant1a, trianglequadDunavant3a, trianglequadDunavant4a, @@ -2090,7 +2090,7 @@ const trianglequadDunavantA = @SVector [ trianglequadDunavant79a ] -const trianglequadDunavantB = @SVector [ +const trianglequadDunavantB = [ trianglequadDunavant1b, trianglequadDunavant3b, trianglequadDunavant4b, @@ -2113,7 +2113,7 @@ const trianglequadDunavantB = @SVector [ trianglequadDunavant79b ] -const trianglequadDunavantW = @SVector [ +const trianglequadDunavantW = [ trianglequadDunavant1w, trianglequadDunavant3w, trianglequadDunavant4w, diff --git a/src/quadrature/triangle/TriangleGauss.jl b/src/quadrature/triangle/TriangleGauss.jl index 94bd6cb..907993b 100644 --- a/src/quadrature/triangle/TriangleGauss.jl +++ b/src/quadrature/triangle/TriangleGauss.jl @@ -5,75 +5,75 @@ # TODO: Add reference to higher order cubature rules -const trianglequadGauss1a = @SVector [ +const trianglequadGauss1a = [ 0.333333333333333, ] -const trianglequadGauss1b = @SVector [ +const trianglequadGauss1b = [ 0.333333333333333, ] -const trianglequadGauss1c = @SVector [ +const trianglequadGauss1c = [ 0.333333333333333, ] -const trianglequadGauss1w = @SVector [ +const trianglequadGauss1w = [ 1.0, ] -const trianglequadGauss3a = @SVector [ +const trianglequadGauss3a = [ 0.666666666666667, 0.166666666666667, 0.166666666666667, ] -const trianglequadGauss3b = @SVector [ +const trianglequadGauss3b = [ 0.166666666666667, 0.666666666666667, 0.166666666666667, ] -const trianglequadGauss3c = @SVector [ +const trianglequadGauss3c = [ 0.166666666666667, 0.166666666666667, 0.666666666666667, ] -const trianglequadGauss3w = @SVector [ +const trianglequadGauss3w = [ 0.333333333333333, 0.333333333333333, 0.333333333333333, ] -const trianglequadGauss4a = @SVector [ +const trianglequadGauss4a = [ 0.333333333333333, 0.6, 0.2, 0.2, ] -const trianglequadGauss4b = @SVector [ +const trianglequadGauss4b = [ 0.333333333333333, 0.2, 0.6, 0.2, ] -const trianglequadGauss4c = @SVector [ +const trianglequadGauss4c = [ 0.333333333333333, 0.2, 0.2, 0.6, ] -const trianglequadGauss4w = @SVector [ +const trianglequadGauss4w = [ -0.562500000000000, 0.520833333333333, 0.520833333333333, 0.520833333333333, ] -const trianglequadGauss6a = @SVector [ +const trianglequadGauss6a = [ 0.816847572980459, 0.091576213509771, 0.091576213509771, @@ -82,7 +82,7 @@ const trianglequadGauss6a = @SVector [ 0.445948490915965, ] -const trianglequadGauss6b = @SVector [ +const trianglequadGauss6b = [ 0.091576213509771, 0.816847572980459, 0.091576213509771, @@ -91,7 +91,7 @@ const trianglequadGauss6b = @SVector [ 0.445948490915965, ] -const trianglequadGauss6c = @SVector [ +const trianglequadGauss6c = [ 0.091576213509771, 0.091576213509771, 0.816847572980459, @@ -100,7 +100,7 @@ const trianglequadGauss6c = @SVector [ 0.108103018168070, ] -const trianglequadGauss6w = @SVector [ +const trianglequadGauss6w = [ 0.109951743655322, 0.109951743655322, 0.109951743655322, @@ -109,7 +109,7 @@ const trianglequadGauss6w = @SVector [ 0.223381589678011, ] -const trianglequadGauss7a = @SVector [ +const trianglequadGauss7a = [ 0.333333333333333, 0.797426985353087, 0.101286507323456, @@ -119,7 +119,7 @@ const trianglequadGauss7a = @SVector [ 0.470142064105115, ] -const trianglequadGauss7b = @SVector [ +const trianglequadGauss7b = [ 0.333333333333333, 0.101286507323456, 0.797426985353087, @@ -129,7 +129,7 @@ const trianglequadGauss7b = @SVector [ 0.470142064105115, ] -const trianglequadGauss7c = @SVector [ +const trianglequadGauss7c = [ 0.333333333333333, 0.101286507323456, 0.101286507323456, @@ -139,7 +139,7 @@ const trianglequadGauss7c = @SVector [ 0.059715871789770, ] -const trianglequadGauss7w = @SVector [ +const trianglequadGauss7w = [ 0.225000000000000, 0.125939180544827, 0.125939180544827, @@ -149,7 +149,7 @@ const trianglequadGauss7w = @SVector [ 0.132394152788506, ] -const trianglequadGauss12a = @SVector [ +const trianglequadGauss12a = [ 0.873821971016996, 0.063089014491502, 0.063089014491502, @@ -164,7 +164,7 @@ const trianglequadGauss12a = @SVector [ 0.053145049844816, ] -const trianglequadGauss12b = @SVector [ +const trianglequadGauss12b = [ 0.063089014491502, 0.873821971016996, 0.063089014491502, @@ -179,7 +179,7 @@ const trianglequadGauss12b = @SVector [ 0.310352451033785, ] -const trianglequadGauss12c = @SVector [ +const trianglequadGauss12c = [ 0.063089014491502, 0.063089014491502, 0.873821971016996, @@ -194,7 +194,7 @@ const trianglequadGauss12c = @SVector [ 0.636502499121399, ] -const trianglequadGauss12w = @SVector [ +const trianglequadGauss12w = [ 0.050844906370207, 0.050844906370207, 0.050844906370207, @@ -209,7 +209,7 @@ const trianglequadGauss12w = @SVector [ 0.082851075618374, ] -const trianglequadGauss13a = @SVector [ +const trianglequadGauss13a = [ 0.333333333333333, 0.479308067841923, 0.260345966079038, @@ -225,7 +225,7 @@ const trianglequadGauss13a = @SVector [ 0.0486903154253160, ] -const trianglequadGauss13b = @SVector [ +const trianglequadGauss13b = [ 0.333333333333333, 0.260345966079038, 0.479308067841923, @@ -241,7 +241,7 @@ const trianglequadGauss13b = @SVector [ 0.312865496004875, ] -const trianglequadGauss13c = @SVector [ +const trianglequadGauss13c = [ 0.333333333333333, 0.260345966079038, 0.260345966079038, @@ -257,7 +257,7 @@ const trianglequadGauss13c = @SVector [ 0.638444188569809, ] -const trianglequadGauss13w = @SVector [ +const trianglequadGauss13w = [ -0.149570044467670, 0.175615257433204, 0.175615257433204, @@ -273,7 +273,7 @@ const trianglequadGauss13w = @SVector [ 0.077113760890257, ] -const trianglequadGauss36a = @SVector [ +const trianglequadGauss36a = [ 0.0242935351590, 0.0265193427722, 0.9492126023551, @@ -312,7 +312,7 @@ const trianglequadGauss36a = @SVector [ 0.2305424298836 ] -const trianglequadGauss36b = @SVector [ +const trianglequadGauss36b = [ 0.9493059293846, 0.0242695130640, 0.0265067966437, @@ -351,7 +351,7 @@ const trianglequadGauss36b = @SVector [ 0.3456013949376 ] -const trianglequadGauss36w = (@SVector [ +const trianglequadGauss36w = ( [ 0.0166240998757, 0.0166811699778, 0.0166830569067, @@ -390,7 +390,7 @@ const trianglequadGauss36w = (@SVector [ 0.1035854367193 ]) ./ 2 -const trianglequadGauss78a = @SVector [ +const trianglequadGauss78a = [ 0.0089411337112, 0.9792622629807, 0.0105475382112, @@ -471,7 +471,7 @@ const trianglequadGauss78a = @SVector [ 0.3729077987144 ] -const trianglequadGauss78b = @SVector [ +const trianglequadGauss78b = [ 0.0086983293702, 0.0102644133744, 0.9785514202515, @@ -552,7 +552,7 @@ const trianglequadGauss78b = @SVector [ 0.3753750277549 ] -const trianglequadGauss78w = (@SVector [ +const trianglequadGauss78w = ( [ 0.0021744545399, 0.0028987135265, 0.0030846029337, @@ -633,7 +633,7 @@ const trianglequadGauss78w = (@SVector [ 0.0637206605672 ]) ./ 2 -const trianglequadGauss105a = @SVector [ +const trianglequadGauss105a = [ 0.0087809303836, 0.9903675314220, 0.0027029276450, @@ -741,7 +741,7 @@ const trianglequadGauss105a = @SVector [ 0.4098894602340 ] -const trianglequadGauss105b = @SVector [ +const trianglequadGauss105b = [ 0.9903676436772, 0.0087809216232, 0.0335914404439, @@ -849,7 +849,7 @@ const trianglequadGauss105b = @SVector [ 0.4098894317792, ] -const trianglequadGauss105w = (@SVector [ +const trianglequadGauss105w = ( [ 0.0006438298261, 0.0006438413076, 0.0010134735710, @@ -957,7 +957,7 @@ const trianglequadGauss105w = (@SVector [ 0.0398766879831 ]) ./ 2 -const trianglequadGauss120a = @SVector [ +const trianglequadGauss120a = [ 0.0082881595033, 0.4618422030241, 0.0071066441239, @@ -1080,7 +1080,7 @@ const trianglequadGauss120a = @SVector [ 0.3361523347440 ] -const trianglequadGauss120b = @SVector [ +const trianglequadGauss120b = [ 0.9848202768869, 0.5381577969759, 0.0080842361390, @@ -1203,7 +1203,7 @@ const trianglequadGauss120b = @SVector [ 0.2778500044356 ] -const trianglequadGauss120w = (@SVector [ +const trianglequadGauss120w = ( [ 0.0014873417859, 0.0014889035262, 0.0015005944380, @@ -1326,7 +1326,7 @@ const trianglequadGauss120w = (@SVector [ 0.0371924811018 ]) ./ 2 -const trianglequadGauss400a = @SVector [ +const trianglequadGauss400a = [ 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, 3.4357004074525577e-03, @@ -1408,7 +1408,7 @@ const trianglequadGauss400a = @SVector [ 9.9656429959254744e-01, 9.9656429959254744e-01, 9.9656429959254744e-01, 9.9656429959254744e-01, 9.9656429959254744e-01, 9.9656429959254744e-01, 9.9656429959254744e-01, 9.9656429959254744e-01, 9.9656429959254744e-01, 9.9656429959254744e-01 ] -const trianglequadGauss400b = @SVector [ +const trianglequadGauss400b = [ 3.4238963701627881e-03, 1.7952145528977594e-02, 4.3732017769028415e-02, 8.0165141146159247e-02, 1.2639828298375833e-01, 1.8134795437803311e-01, 2.4372624180423841e-01, 3.1207107651919747e-01, 3.8478052619624081e-01, 4.6015035032944468e-01, 5.3641394926310271e-01, 6.1178377339630663e-01, 6.8449322307334992e-01, 7.5283805778830903e-01, 8.1521634521451436e-01, @@ -1490,7 +1490,7 @@ const trianglequadGauss400b = @SVector [ 1.8493113036459207e-03, 2.1091521745159178e-03, 2.3598212843598133e-03, 2.5954431871045409e-03, 2.8104951487431724e-03, 2.9999366212862710e-03, 3.1593274647212558e-03, 3.2849323021439448e-03, 3.3738095753870568e-03, 3.4238963701627881e-03 ] -const trianglequadGauss400w = (@SVector [ +const trianglequadGauss400w = ( [ 3.0918731028920757e-04, 7.1269681990617786e-04, 1.1001132168084852e-03, 1.4617975077641802e-03, 1.7892294090025515e-03, 2.0747266161171671e-03, 2.3115952886863795e-03, 2.4942827316709218e-03, 2.6185066288951639e-03, 2.6813551585055523e-03, 2.6813551585055523e-03, 2.6185066288951639e-03, 2.4942827316709218e-03, 2.3115952886863795e-03, 2.0747266161171671e-03, @@ -1572,7 +1572,7 @@ const trianglequadGauss400w = (@SVector [ 9.2440929444984856e-06, 9.0274198017032943e-06, 8.5991523085940462e-06, 7.9693291024496169e-06, 7.1527136615879260e-06, 6.1684491528038026e-06, 5.0396128931087845e-06, 3.7926899737209111e-06, 2.4570544575427101e-06, 1.0659372088425298e-06 ]) ./ 2 -const trianglequadGauss900a = @SVector [ +const trianglequadGauss900a = [ 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, 1.5532579626752474e-03, @@ -1754,7 +1754,7 @@ const trianglequadGauss900a = @SVector [ 9.9844674203732475e-01, 9.9844674203732475e-01, 9.9844674203732475e-01, 9.9844674203732475e-01, 9.9844674203732475e-01, 9.9844674203732475e-01, 9.9844674203732475e-01, 9.9844674203732475e-01, 9.9844674203732475e-01, 9.9844674203732475e-01 ] -const trianglequadGauss900b = @SVector [ +const trianglequadGauss900b = [ 1.5508453523766334e-03, 8.1532545513458305e-03, 1.9958019337560785e-02, 3.6842661103375078e-02, 5.8628525212516977e-02, 8.5084754640270099e-02, 1.1593093317123887e-01, 1.5084009454090175e-01, 1.8944219814142918e-01, 2.3132805481318516e-01, 2.7605366563067762e-01, 3.2314492896492530e-01, 3.7210266635129818e-01, 4.2240791404863498e-01, 4.7352742426565486e-01, @@ -1936,7 +1936,7 @@ const trianglequadGauss900b = @SVector [ 1.1933868468703306e-03, 1.2585475987258832e-03, 1.3185999002349133e-03, 1.3729071863272115e-03, 1.4208937943295222e-03, 1.4620510712185931e-03, 1.4959427806874798e-03, 1.5222097843898060e-03, 1.5405741538946657e-03, 1.5508453523766334e-03 ] -const trianglequadGauss900w = (@SVector [ +const trianglequadGauss900w = ( [ 6.3393472058975955e-05, 1.4691582105612550e-04, 2.2900583486229500e-04, 3.0867923073618977e-04, 3.8508275109363019e-04, 4.5740496279664960e-04, 5.2487882297804151e-04, 5.8678895194472442e-04, 6.4247903025549066e-04, 6.9135870330888595e-04, 7.3290982143010656e-04, 7.6669192547512504e-04, 7.9234691283714938e-04, 8.0960283204461417e-04, 8.1827676484210919e-04, @@ -2118,7 +2118,7 @@ const trianglequadGauss900w = (@SVector [ 1.0755289849393060e-06, 9.9948813249661188e-07, 9.1285250744394168e-07, 8.1654050927813780e-07, 7.1157315730367923e-07, 5.9906334934258464e-07, 4.8020435429046972e-07, 3.5625849784750484e-07, 2.2855317093104202e-07, 9.8619596931444806e-08 ]) ./ 2 -const trianglequadGaussA = @SVector [ +const trianglequadGaussA = [ trianglequadGauss1a, trianglequadGauss3a, trianglequadGauss4a, @@ -2134,7 +2134,7 @@ const trianglequadGaussA = @SVector [ trianglequadGauss900a ] -const trianglequadGaussB = @SVector [ +const trianglequadGaussB = [ trianglequadGauss1b, trianglequadGauss3b, trianglequadGauss4b, @@ -2150,7 +2150,7 @@ const trianglequadGaussB = @SVector [ trianglequadGauss900b ] -const trianglequadGaussW = @SVector [ +const trianglequadGaussW = [ trianglequadGauss1w, trianglequadGauss3w, trianglequadGauss4w, diff --git a/src/quadrature/triangle/TriangleQuadratures.jl b/src/quadrature/triangle/TriangleQuadratures.jl index b6f564d..d549dcf 100644 --- a/src/quadrature/triangle/TriangleQuadratures.jl +++ b/src/quadrature/triangle/TriangleQuadratures.jl @@ -27,9 +27,19 @@ function trgauss(n) end function trgauss(trqd::TriangleQuadDunavant) + # n = trqd.order + # @assert 1 <= n <= length(trianglequadDunavantA) + # u = copy(transpose([trianglequadDunavantA[n] trianglequadDunavantB[n]])) + # w = trianglequadDunavantW[n] + # return u, w/2 + n = trqd.order @assert 1 <= n <= length(trianglequadDunavantA) - u = copy(transpose([trianglequadDunavantA[n] trianglequadDunavantB[n]])) + A = trianglequadDunavantA[n] + B = trianglequadDunavantB[n] + T = eltype(A) + AB = [A, B] + u = T[AB[i][j] for i in (1,2), j in 1:length(A)] w = trianglequadDunavantW[n] return u, w/2 end @@ -37,7 +47,11 @@ end function trgauss(trqd::TriangleQuadLegacy) n = trqd.order @assert 1 <= n <= length(trianglequadGaussA) - u = copy(transpose([trianglequadGaussA[n] trianglequadGaussB[n]])) + A = trianglequadGaussA[n] + B = trianglequadGaussB[n] + T = eltype(A) + AB = [A, B] + u = T[AB[i][j] for i in (1,2), j in 1:length(A)] w = trianglequadGaussW[n] return u, w/2 end \ No newline at end of file