From 732595d46bd7f0c37e91da31c4469fbce2f09a00 Mon Sep 17 00:00:00 2001 From: Anthony Royer <anthony.royer@uliege.be> Date: Thu, 10 Nov 2022 15:36:39 +0100 Subject: [PATCH] New model acoustic --- .../helmholtz/crossPoints/Subproblem3D.cpp | 92 +++++++++++++- examples/helmholtz/crossPoints/Subproblem3D.h | 18 +++ examples/helmholtz/crossPoints/ddm3D.cpp | 107 +++++++++++----- examples/helmholtz/crossPoints/mesh.cpp | 2 + examples/navier/crossPoints/Subproblem2D.cpp | 2 +- examples/navier/crossPoints/Subproblem3D.cpp | 117 +++++++++--------- examples/navier/crossPoints/ddm3D.cpp | 8 +- examples/navier/crossPoints/mesh.cpp | 5 +- 8 files changed, 248 insertions(+), 103 deletions(-) diff --git a/examples/helmholtz/crossPoints/Subproblem3D.cpp b/examples/helmholtz/crossPoints/Subproblem3D.cpp index 0b1f042e..0517c841 100755 --- a/examples/helmholtz/crossPoints/Subproblem3D.cpp +++ b/examples/helmholtz/crossPoints/Subproblem3D.cpp @@ -240,7 +240,8 @@ namespace D3 { _corner[c] = new HABC_HABC_HABC(c, static_cast< HABC_HABC * >(_edge[c%4]), static_cast< HABC_HABC * >(_edge[(c+1)%4]), static_cast< HABC_HABC * >(_edge[c+4])); } else if(method[c%4] == "pmlContinuous" && method[(c+1)%4] == "pmlContinuous" && method[4] == "pmlContinuous") { - _corner[c] = new PmlContinuous_PmlContinuous_PmlContinuous(c, static_cast< PmlContinuous_PmlContinuous * >(_edge[c%4]), static_cast< PmlContinuous_PmlContinuous * >(_edge[(c+1)%4]), static_cast< PmlContinuous_PmlContinuous * >(_edge[c+4])); + _corner[c] = new PmlContinuous_PmlContinuous_PmlContinuous_simplify(c, static_cast< PmlContinuous_PmlContinuous * >(_edge[c%4]), static_cast< PmlContinuous_PmlContinuous * >(_edge[(c+1)%4]), static_cast< PmlContinuous_PmlContinuous * >(_edge[c+4])); +// _corner[c] = new PmlContinuous_PmlContinuous_PmlContinuous(c, static_cast< PmlContinuous_PmlContinuous * >(_edge[c%4]), static_cast< PmlContinuous_PmlContinuous * >(_edge[(c+1)%4]), static_cast< PmlContinuous_PmlContinuous * >(_edge[c+4])); } } else { @@ -248,7 +249,8 @@ namespace D3 { _corner[c] = new HABC_HABC_HABC(c, static_cast< HABC_HABC * >(_edge[c%4+8]), static_cast< HABC_HABC * >(_edge[(c+1)%4+8]), static_cast< HABC_HABC * >(_edge[c])); } else if(method[c%4] == "pmlContinuous" && method[(c+1)%4] == "pmlContinuous" && method[5] == "pmlContinuous") { - _corner[c] = new PmlContinuous_PmlContinuous_PmlContinuous(c, static_cast< PmlContinuous_PmlContinuous * >(_edge[c%4+8]), static_cast< PmlContinuous_PmlContinuous * >(_edge[(c+1)%4+8]), static_cast< PmlContinuous_PmlContinuous * >(_edge[c])); + _corner[c] = new PmlContinuous_PmlContinuous_PmlContinuous_simplify(c, static_cast< PmlContinuous_PmlContinuous * >(_edge[c%4+8]), static_cast< PmlContinuous_PmlContinuous * >(_edge[(c+1)%4+8]), static_cast< PmlContinuous_PmlContinuous * >(_edge[c])); +// _corner[c] = new PmlContinuous_PmlContinuous_PmlContinuous(c, static_cast< PmlContinuous_PmlContinuous * >(_edge[c%4+8]), static_cast< PmlContinuous_PmlContinuous * >(_edge[(c+1)%4+8]), static_cast< PmlContinuous_PmlContinuous * >(_edge[c])); } } } @@ -491,7 +493,7 @@ namespace D3 { sigma = sigmaStar * distSigma * distSigma / (_size * _size); } - gmshfem::function::ScalarFunction< std::complex< double > > kMod = 1. + im * sigma / kappaMod; + gmshfem::function::ScalarFunction< std::complex< double > > kMod = 1. + im * sigma/* / kappaMod*/; if(orientation() == "E" || orientation() == "W") { D = gmshfem::function::tensorDiag< std::complex< double > >(1./kMod, kMod, kMod); E = kMod; @@ -553,10 +555,10 @@ namespace D3 { formulation.integral(- D * grad(dof(*_uPml)), grad(tf(*_uPml)), pml, gauss); formulation.integral(kappaMod * kappaMod * E * dof(*_uPml), tf(*_uPml), pml, gauss); - + formulation.integral(dof(*_v), tf(*_uPml), sigma, gauss); formulation.integral(- dof(*_v), tf(*u), sigma, gauss); - + formulation.integral(dof(*_uPml), tf(*_v), sigma, gauss); formulation.integral(- dof(*u), tf(*_v), sigma, gauss); } @@ -1331,6 +1333,86 @@ namespace D3 { formulation.integral(dof(*_vEdge[2]), tf(*_vCorner), corner, gauss); formulation.integral(dof(*_vCorner), tf(*_vEdge[2]), corner, gauss); } + + // ********************************** + // PML continuous _ PML continuous _ PML continuous simplify + // ********************************** + + PmlContinuous_PmlContinuous_PmlContinuous_simplify::PmlContinuous_PmlContinuous_PmlContinuous_simplify(const unsigned int corner, PmlContinuous_PmlContinuous *first, PmlContinuous_PmlContinuous *second, PmlContinuous_PmlContinuous *third) : Corner(corner, first, second, third) + { + } + + PmlContinuous_PmlContinuous_PmlContinuous_simplify::~PmlContinuous_PmlContinuous_PmlContinuous_simplify() + { + delete _vC[0]; + delete _vC[1]; + delete _vC[2]; + } + + gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* PmlContinuous_PmlContinuous_PmlContinuous_simplify::getV(const unsigned int i) const + { + return _vC[i]; + } + + PmlContinuous_PmlContinuous *PmlContinuous_PmlContinuous_PmlContinuous_simplify::firstEdge() const + { + return static_cast< PmlContinuous_PmlContinuous * >(_bnd[0]); + } + + PmlContinuous_PmlContinuous *PmlContinuous_PmlContinuous_PmlContinuous_simplify::secondEdge() const + { + return static_cast< PmlContinuous_PmlContinuous * >(_bnd[1]); + } + + PmlContinuous_PmlContinuous *PmlContinuous_PmlContinuous_PmlContinuous_simplify::thirdEdge() const + { + return static_cast< PmlContinuous_PmlContinuous * >(_bnd[2]); + } + + void PmlContinuous_PmlContinuous_PmlContinuous_simplify::writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters ¶meters) + { + gmshfem::domain::Domain corner = domains.getCorner(_corner); + gmshfem::domain::Domain pmlCorner = domains.getPmlCorner(_corner); + std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > pmlBnd = domains.getPmlEdgeBnd(_corner); + std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > cornerEdge = domains.getCornerEdge(_corner); + + const unsigned int neumannOrder = parameters.getNeumannOrder(); + const unsigned int fieldOrder = parameters.getFieldOrder(); + const std::string gauss = parameters.getGauss(); + const gmshfem::function::ScalarFunction< std::complex< double > > kappa = parameters.getKappa(); + + _vC[0] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vC_" + orientation() + "_first", std::get<0>(pmlBnd), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); + _vC[1] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vC_" + orientation() + "_second", std::get<1>(pmlBnd), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); + _vC[2] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vC_" + orientation() + "_third", std::get<2>(pmlBnd), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); + + gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *uPml[3]; + uPml[0] = static_cast< PmlContinuous_PmlContinuous * >(_bnd[0])->getUPml(); + uPml[1] = static_cast< PmlContinuous_PmlContinuous * >(_bnd[1])->getUPml(); + uPml[2] = static_cast< PmlContinuous_PmlContinuous * >(_bnd[2])->getUPml(); + + std::vector< PmlContinuous * > boundaries(6); + for(unsigned int i = 0; i < 3; ++i) { + boundaries[2 * i] = static_cast< PmlContinuous_PmlContinuous * >(_bnd[i])->firstBoundary(); + boundaries[2 * i + 1] = static_cast< PmlContinuous_PmlContinuous * >(_bnd[i])->secondBoundary(); + } + std::sort(boundaries.begin(), boundaries.end()); + gmshfem::function::TensorFunction< std::complex< double > > D[3]; + gmshfem::function::ScalarFunction< std::complex< double > > E[3]; + gmshfem::function::ScalarFunction< std::complex< double > > kappaMod = boundaries[0]->pmlCoefficients(D[0], E[0], corner, kappa); + boundaries[2]->pmlCoefficients(D[1], E[1], corner, kappa); + boundaries[4]->pmlCoefficients(D[2], E[2], corner, kappa); + gmshfem::function::ScalarFunction< std::complex< double > > Et = E[0] * E[1] * E[2]; + + gmshfem::domain::Domain pmlBnds[3] {std::get<0>(pmlBnd), std::get<1>(pmlBnd), std::get<2>(pmlBnd)}; + const std::complex< double > im(0., 1.); + + for(unsigned int i = 0; i < 3; ++i) { + formulation.integral(- dof(*_vC[i]), tf(*uPml[i]), pmlBnds[i], gauss); + + formulation.integral(dof(*_vC[i]), tf(*_vC[i]), pmlBnds[i], gauss); + formulation.integral(im * kappaMod * Et * dof(*uPml[i]), tf(*_vC[i]), pmlBnds[i], gauss); + } + } } diff --git a/examples/helmholtz/crossPoints/Subproblem3D.h b/examples/helmholtz/crossPoints/Subproblem3D.h index c1c27170..5b3fd281 100755 --- a/examples/helmholtz/crossPoints/Subproblem3D.h +++ b/examples/helmholtz/crossPoints/Subproblem3D.h @@ -266,6 +266,24 @@ namespace D3 { void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters ¶meters) override; }; + + class PmlContinuous_PmlContinuous_PmlContinuous_simplify : public Corner + { + protected: + gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _vC[3]; + + public: + PmlContinuous_PmlContinuous_PmlContinuous_simplify(const unsigned int corner, PmlContinuous_PmlContinuous *first, PmlContinuous_PmlContinuous *second, PmlContinuous_PmlContinuous *third); + ~PmlContinuous_PmlContinuous_PmlContinuous_simplify(); + + gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getV(const unsigned int i) const; + + PmlContinuous_PmlContinuous *firstEdge() const; + PmlContinuous_PmlContinuous *secondEdge() const; + PmlContinuous_PmlContinuous *thirdEdge() const; + + void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters ¶meters) override; + }; } diff --git a/examples/helmholtz/crossPoints/ddm3D.cpp b/examples/helmholtz/crossPoints/ddm3D.cpp index d82e5ecb..60eb8bf7 100755 --- a/examples/helmholtz/crossPoints/ddm3D.cpp +++ b/examples/helmholtz/crossPoints/ddm3D.cpp @@ -41,7 +41,7 @@ namespace D3 { // P H Y S I C S // ************************ double pi = 3.14159265359; - double k = 0.5 * pi; + double k = 1. * pi; gmshDdm->userDefinedParameter(k, "k"); double R = 0.5; gmshDdm->userDefinedParameter(R, "R"); @@ -51,7 +51,7 @@ namespace D3 { // ************************ // M E S H // ************************ - double lc = 0.4; // other = 0.06666666666666 + double lc = 0.2; // other = 0.06666666666666 gmshDdm->userDefinedParameter(lc, "lc"); int meshOrder = 1; gmshDdm->userDefinedParameter(meshOrder, "meshOrder"); @@ -467,22 +467,19 @@ namespace D3 { } } + std::vector< int > sourcesDom; + for(unsigned int i = 0; i < nDomX-1; ++i) { + for(unsigned int j = 0; j < nDomY-1; ++j) { + sourcesDom.push_back(i * nDomY * nDomZ + j * nDomZ + nDomZ-1); + } + } if(benchmark == "scattering") { gamma = gmshfem::domain::Domain("gamma"); } else if(benchmark == "salt") { - gamma = corner[4](1) | - corner[4](3) | - corner[4](5) | - corner[4](9) | - corner[4](11) | - corner[4](13) | - corner[4](17) | - corner[4](19) | - corner[4](21) | - corner[4](25) | - corner[4](27) | - corner[4](29); + for(auto i : sourcesDom) { + gamma |= corner[4](i); + } } // Kappa definition @@ -537,7 +534,7 @@ namespace D3 { const double max = *std::max_element(dataRaw.begin(), dataRaw.end()); gmshfem::msg::info << "min = " << min << ", max = " << max << gmshfem::msg::endl; - kappa = 2. * pi * 1.3 / gmshfem::function::trilinearInterpolation< std::complex< double > >(&x, &y, &z, &data); + kappa = 2. * pi * 4. / gmshfem::function::trilinearInterpolation< std::complex< double > >(&x, &y, &z, &data); for(unsigned int k = 0; k < nDomX * nDomY * nDomZ; ++k) { if(gmshddm::mpi::isItMySubdomain(k)) { //gmshfem::post::save(kappa, omega(k), "k_" + std::to_string(k), "msh"); @@ -551,18 +548,9 @@ namespace D3 { u(0).domain(u(0).domain() | gamma); } else if (benchmark == "salt") { - u(1).addConstraint(corner[4](1), 1.); - u(3).addConstraint(corner[4](3), 1.); - u(5).addConstraint(corner[4](5), 1.); - u(9).addConstraint(corner[4](9), 1.); - u(11).addConstraint(corner[4](11), 1.); - u(13).addConstraint(corner[4](13), 1.); - u(17).addConstraint(corner[4](17), 1.); - u(19).addConstraint(corner[4](19), 1.); - u(21).addConstraint(corner[4](21), 1.); - u(25).addConstraint(corner[4](25), 1.); - u(27).addConstraint(corner[4](27), 1.); - u(29).addConstraint(corner[4](29), 1.); + for(auto i : sourcesDom) { + u(i).addConstraint(corner[4](i), 1.); + } } gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > lambda("lambda", gamma, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder); @@ -771,7 +759,7 @@ namespace D3 { if(index == 0) { formulation(index).integral(dof(lambda), tf(u(index)), gamma, gauss); formulation(index).integral(dof(u(index)), tf(lambda), gamma, gauss); - formulation(index).integral(formulation.physicalSource(- fAnalytic), tf(lambda), gamma, gauss); + formulation(index).integral(formulation.physicalSource(- 1.), tf(lambda), gamma, gauss); } } @@ -960,6 +948,20 @@ namespace D3 { formulation.artificialSourceTerm(termId); } } + else if(auto cr = dynamic_cast< PmlContinuous_PmlContinuous_PmlContinuous_simplify * >(corner)) { + if(!std::get<0>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) { + auto termId = formulation(index).integral(- (*std::get<1>(gCornerPmlContinuous[cornerNeighbor[c][0]]))(jndex, index), tf(*cr->firstEdge()->getUPml()), std::get<0>(pmlEdgeBndInterface[c])(index, jndex), gauss); + formulation.artificialSourceTerm(termId); + } + if(!std::get<1>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) { + auto termId = formulation(index).integral(- (*std::get<0>(gCornerPmlContinuous[cornerNeighbor[c][1]]))(jndex, index), tf(*cr->secondEdge()->getUPml()), std::get<1>(pmlEdgeBndInterface[c])(index, jndex), gauss); + formulation.artificialSourceTerm(termId); + } + if(!std::get<2>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) { + auto termId = formulation(index).integral(- (*std::get<2>(gCornerPmlContinuous[cornerNeighbor[c][2]]))(jndex, index), tf(*cr->thirdEdge()->getUPml()), std::get<2>(pmlEdgeBndInterface[c])(index, jndex), gauss); + formulation.artificialSourceTerm(termId); + } + } } } @@ -1135,6 +1137,26 @@ namespace D3 { formulation(index, jndex).integral(2. * *cr->getV(2), tf((*std::get<2>(gCornerPmlContinuous[c]))(index, jndex)), std::get<2>(pmlEdgeBndInterface[c])(index, jndex), gauss); } } + else if(auto cr = dynamic_cast< PmlContinuous_PmlContinuous_PmlContinuous_simplify * >(corner)) { + if(!std::get<0>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) { + formulation(index, jndex).integral(dof((*std::get<0>(gCornerPmlContinuous[c]))(index, jndex)), tf((*std::get<0>(gCornerPmlContinuous[c]))(index, jndex)), std::get<0>(pmlEdgeBndInterface[c])(index, jndex), gauss); + auto termId = formulation(index, jndex).integral((*std::get<1>(gCornerPmlContinuous[cornerNeighbor[c][0]]))(jndex, index), tf((*std::get<0>(gCornerPmlContinuous[c]))(index, jndex)), std::get<0>(pmlEdgeBndInterface[c])(index, jndex), gauss); + formulation.artificialSourceTerm(termId); + formulation(index, jndex).integral(2. * *cr->getV(0), tf((*std::get<0>(gCornerPmlContinuous[c]))(index, jndex)), std::get<0>(pmlEdgeBndInterface[c])(index, jndex), gauss); + } + if(!std::get<1>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) { + formulation(index, jndex).integral(dof((*std::get<1>(gCornerPmlContinuous[c]))(index, jndex)), tf((*std::get<1>(gCornerPmlContinuous[c]))(index, jndex)), std::get<1>(pmlEdgeBndInterface[c])(index, jndex), gauss); + auto termId = formulation(index, jndex).integral((*std::get<0>(gCornerPmlContinuous[cornerNeighbor[c][1]]))(jndex, index), tf((*std::get<1>(gCornerPmlContinuous[c]))(index, jndex)), std::get<1>(pmlEdgeBndInterface[c])(index, jndex), gauss); + formulation.artificialSourceTerm(termId); + formulation(index, jndex).integral(2. * *cr->getV(1), tf((*std::get<1>(gCornerPmlContinuous[c]))(index, jndex)), std::get<1>(pmlEdgeBndInterface[c])(index, jndex), gauss); + } + if(!std::get<2>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) { + formulation(index, jndex).integral(dof((*std::get<2>(gCornerPmlContinuous[c]))(index, jndex)), tf((*std::get<2>(gCornerPmlContinuous[c]))(index, jndex)), std::get<2>(pmlEdgeBndInterface[c])(index, jndex), gauss); + auto termId = formulation(index, jndex).integral((*std::get<2>(gCornerPmlContinuous[cornerNeighbor[c][2]]))(jndex, index), tf((*std::get<2>(gCornerPmlContinuous[c]))(index, jndex)), std::get<2>(pmlEdgeBndInterface[c])(index, jndex), gauss); + formulation.artificialSourceTerm(termId); + formulation(index, jndex).integral(2. * *cr->getV(2), tf((*std::get<2>(gCornerPmlContinuous[c]))(index, jndex)), std::get<2>(pmlEdgeBndInterface[c])(index, jndex), gauss); + } + } } } } @@ -1142,8 +1164,8 @@ namespace D3 { } formulation.pre(); - formulation.solve("gmres", res, iterMax); + gmshfem::msg::error << "coucou" << gmshfem::msg::endl; for(unsigned int i = 0; i < nDomX; ++i) { for(unsigned int j = 0; j < nDomY; ++j) { for(unsigned int k = 0; k < nDomZ; ++k) { @@ -1156,6 +1178,33 @@ namespace D3 { } } } + + formulation.solve("gmres", res, iterMax); + +// { +// Boundary *boundary = subproblem[0][0][0]->getBoundary(0); +// auto bd = dynamic_cast< PmlContinuous * >(boundary); +// gmshfem::post::save( *bd->getV(), bd->getV()->domain(), bd->getV()->name()); +// } +// +// { +// Corner *corner = subproblem[0][0][0]->getCorner(4); +// auto cr = dynamic_cast< PmlContinuous_PmlContinuous_PmlContinuous_simplify * >(corner); +// gmshfem::post::save( *cr->getV(2), cr->getV(2)->domain(), cr->getV(2)->name()); +// } + +// for(unsigned int i = 0; i < nDomX; ++i) { +// for(unsigned int j = 0; j < nDomY; ++j) { +// for(unsigned int k = 0; k < nDomZ; ++k) { +// unsigned int index = i * nDomY * nDomZ + j * nDomZ + k; +// if(gmshddm::mpi::isItMySubdomain(index)) { +// if(saveU) { +// gmshfem::post::save(u(index), omega(index), "u_" + std::to_string(index), "msh"); +// } +// } +// } +// } +// } if(monoDomainError) { gmshfem::problem::Formulation< std::complex< double > > formulationMono("HelmholtzMonoDomain"); diff --git a/examples/helmholtz/crossPoints/mesh.cpp b/examples/helmholtz/crossPoints/mesh.cpp index 37197a26..746c859e 100755 --- a/examples/helmholtz/crossPoints/mesh.cpp +++ b/examples/helmholtz/crossPoints/mesh.cpp @@ -1362,6 +1362,7 @@ namespace D3 { void checkerboard(const unsigned int nDomX, const unsigned int nDomY, const unsigned int nDomZ, const double xSize, const double ySize, const double zSize, const double sphereSize, const double lc, const bool pml, const double pmlSize, const bool pmlExt, const double pmlSizeExt, const int order, const unsigned int spherePosX, const unsigned int spherePosY, const unsigned int spherePosZ, bool transfinite) { + int subdomainID = 0; for(int i = 0; i < nDomX; ++i) { for(int j = 0; j < nDomY; ++j) { for(int k = 0; k < nDomZ; ++k) { @@ -1438,6 +1439,7 @@ namespace D3 { } subdomain(i, j, k, xSize, ySize, zSize, (spherePosX == i && spherePosY == j && spherePosZ == k ? sphereSize : 0.), lc, pmlSizeE, pmlSizeN, pmlSizeW, pmlSizeS, pmlSizeD, pmlSizeU, transfinite); + gmshfem::msg::info << "Generate subdomain " << ++subdomainID << gmshfem::msg::endl; } } } diff --git a/examples/navier/crossPoints/Subproblem2D.cpp b/examples/navier/crossPoints/Subproblem2D.cpp index 966d1cbd..11d03d4c 100755 --- a/examples/navier/crossPoints/Subproblem2D.cpp +++ b/examples/navier/crossPoints/Subproblem2D.cpp @@ -1189,7 +1189,7 @@ namespace D2 { formulation.integral((im/kappaS) * It * dof(*_uSymSymPmlCorner), tf(*_uSymSymPmlCorner), pmlCornerInf, gauss); // corner equation - _vSymSymCorner = new gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 2 >("vSymSymCorner_" + orientation() + "_second", corner, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); + _vSymSymCorner = new gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 2 >("vSymSymCorner_" + orientation(), corner, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); formulation.integral(dof(*_vSymSymC[0]), tf(*_vSymSymCorner), corner, gauss); formulation.integral(dof(*_vSymSymCorner), tf(*_vSymSymC[0]), corner, gauss); diff --git a/examples/navier/crossPoints/Subproblem3D.cpp b/examples/navier/crossPoints/Subproblem3D.cpp index 5b3500ef..fa1f03bf 100755 --- a/examples/navier/crossPoints/Subproblem3D.cpp +++ b/examples/navier/crossPoints/Subproblem3D.cpp @@ -234,8 +234,8 @@ namespace D3 { _corner[c] = new PmlContinuous_PmlContinuous_PmlContinuous(c, static_cast< PmlContinuous_PmlContinuous * >(_edge[c%4]), static_cast< PmlContinuous_PmlContinuous * >(_edge[(c+1)%4]), static_cast< PmlContinuous_PmlContinuous * >(_edge[c+4])); } else if(method[c%4] == "pmlContinuousSymmetric" && method[(c+1)%4] == "pmlContinuousSymmetric" && method[4] == "pmlContinuousSymmetric") { - _corner[c] = new PmlContinuousSymmetric_PmlContinuousSymmetric_PmlContinuousSymmetric_simplify(c, static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[c%4]), static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[(c+1)%4]), static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[c+4])); -// _corner[c] = new PmlContinuousSymmetric_PmlContinuousSymmetric_PmlContinuousSymmetric(c, static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[c%4]), static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[(c+1)%4]), static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[c+4])); +// _corner[c] = new PmlContinuousSymmetric_PmlContinuousSymmetric_PmlContinuousSymmetric_simplify(c, static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[c%4]), static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[(c+1)%4]), static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[c+4])); + _corner[c] = new PmlContinuousSymmetric_PmlContinuousSymmetric_PmlContinuousSymmetric(c, static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[c%4]), static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[(c+1)%4]), static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[c+4])); } } else { @@ -243,8 +243,8 @@ namespace D3 { _corner[c] = new PmlContinuous_PmlContinuous_PmlContinuous(c, static_cast< PmlContinuous_PmlContinuous * >(_edge[c%4+8]), static_cast< PmlContinuous_PmlContinuous * >(_edge[(c+1)%4+8]), static_cast< PmlContinuous_PmlContinuous * >(_edge[c])); } else if(method[c%4] == "pmlContinuousSymmetric" && method[(c+1)%4] == "pmlContinuousSymmetric" && method[5] == "pmlContinuousSymmetric") { - _corner[c] = new PmlContinuousSymmetric_PmlContinuousSymmetric_PmlContinuousSymmetric_simplify(c, static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[c%4+8]), static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[(c+1)%4+8]), static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[c])); -// _corner[c] = new PmlContinuousSymmetric_PmlContinuousSymmetric_PmlContinuousSymmetric(c, static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[c%4+8]), static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[(c+1)%4+8]), static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[c])); +// _corner[c] = new PmlContinuousSymmetric_PmlContinuousSymmetric_PmlContinuousSymmetric_simplify(c, static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[c%4+8]), static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[(c+1)%4+8]), static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[c])); + _corner[c] = new PmlContinuousSymmetric_PmlContinuousSymmetric_PmlContinuousSymmetric(c, static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[c%4+8]), static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[(c+1)%4+8]), static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_edge[c])); } } } @@ -563,7 +563,7 @@ namespace D3 { tmpY(i,j)(k,l) *= -1.; // symmetry of grad_y } else if(l == 2) { - tmpY(i,j)(k,l) *= -1.; // symmetry of grad_z + tmpZ(i,j)(k,l) *= -1.; // symmetry of grad_z } } } @@ -975,11 +975,11 @@ namespace D3 { formulation.integral(-dof(*_vE[1]), tf(*_vEdge), edge, gauss); formulation.integral(-dof(*_vEdge), tf(*_vE[1]), edge, gauss); - formulation.integral(dof(*static_cast< PmlContinuous * >(_bnd[0])->getV()), tf(*_vEdge), edge, gauss); - formulation.integral(dof(*_vEdge), tf(*static_cast< PmlContinuous * >(_bnd[0])->getV()), edge, gauss); + formulation.integral(dof(*static_cast< PmlContinuousSymmetric * >(_bnd[0])->getV()), tf(*_vEdge), edge, gauss); + formulation.integral(dof(*_vEdge), tf(*static_cast< PmlContinuousSymmetric * >(_bnd[0])->getV()), edge, gauss); - formulation.integral(-dof(*static_cast< PmlContinuous * >(_bnd[1])->getV()), tf(*_vEdge), edge, gauss); - formulation.integral(-dof(*_vEdge), tf(*static_cast< PmlContinuous * >(_bnd[1])->getV()), edge, gauss); + formulation.integral(-dof(*static_cast< PmlContinuousSymmetric * >(_bnd[1])->getV()), tf(*_vEdge), edge, gauss); + formulation.integral(-dof(*_vEdge), tf(*static_cast< PmlContinuousSymmetric * >(_bnd[1])->getV()), edge, gauss); // symmetric bnd 0 @@ -1112,7 +1112,7 @@ namespace D3 { formulation.integral((im/kappaS) * It * dof(*_uSymSymPmlEdge), tf(*_uSymSymPmlEdge), pmlEdgeInf, gauss); // corner equation - _vSymSymEdge = new gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 >("vSymSymEdge_" + orientation() + "_second", edge, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); + _vSymSymEdge = new gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 >("vSymSymEdge_" + orientation(), edge, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); formulation.integral(dof(*_vSymSymE[0]), tf(*_vSymSymEdge), edge, gauss); formulation.integral(dof(*_vSymSymEdge), tf(*_vSymSymE[0]), edge, gauss); @@ -1291,7 +1291,7 @@ namespace D3 { } // ********************************** - // PML continuous symmetric _ PML continuous symmetric _ PML continuous symmetric + // PML continuous symmetric _ PML continuous symmetric _ PML continuous symmetric _ simplify // ********************************** PmlContinuousSymmetric_PmlContinuousSymmetric_PmlContinuousSymmetric_simplify::PmlContinuousSymmetric_PmlContinuousSymmetric_PmlContinuousSymmetric_simplify(const unsigned int corner, PmlContinuousSymmetric_PmlContinuousSymmetric *first, PmlContinuousSymmetric_PmlContinuousSymmetric *second, PmlContinuousSymmetric_PmlContinuousSymmetric *third) : Corner(corner, first, second, third) @@ -1300,18 +1300,17 @@ namespace D3 { PmlContinuousSymmetric_PmlContinuousSymmetric_PmlContinuousSymmetric_simplify::~PmlContinuousSymmetric_PmlContinuousSymmetric_PmlContinuousSymmetric_simplify() { - delete _vC[0]; - delete _vC[1]; - delete _vC[2]; - delete _vSymC[0][0]; - delete _vSymC[0][1]; - delete _vSymC[0][2]; - delete _vSymC[1][0]; - delete _vSymC[1][1]; - delete _vSymC[1][2]; - delete _vSymSymC[0]; - delete _vSymSymC[1]; - delete _vSymSymC[2]; + for(int i = 0; i < 3; ++i) { + delete _vC[i]; + } + for(int i = 0; i < 3; ++i) { + for(int j = 0; j < 2; ++j) { + delete _vSymC[j][i]; + } + } + for(int i = 0; i < 3; ++i) { + delete _vSymSymC[i]; + } } gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 >* PmlContinuousSymmetric_PmlContinuousSymmetric_PmlContinuousSymmetric_simplify::getV(const unsigned int i) const @@ -1359,42 +1358,40 @@ namespace D3 { gmshfem::function::ScalarFunction< std::complex< double > > kappaS = parameters.getKappaS(); gmshfem::function::ScalarFunction< std::complex< double > > fake = 1.; - _vC[0] = new gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 >("vC_" + orientation() + "_first", std::get<0>(pmlBnd), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); - _vC[1] = new gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 >("vC_" + orientation() + "_second", std::get<1>(pmlBnd), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); - _vC[2] = new gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 >("vC_" + orientation() + "_third", std::get<2>(pmlBnd), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); - - _vSymC[0][0] = new gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 >("vSymC_first_" + orientation() + "_first", std::get<0>(pmlBnd), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); - _vSymC[0][1] = new gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 >("vSymC_first_" + orientation() + "_second", std::get<1>(pmlBnd), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); - _vSymC[0][2] = new gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 >("vSymC_first_" + orientation() + "_third", std::get<2>(pmlBnd), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); + std::string IlearnToCount[3] {"first", "second", "third"}; + gmshfem::domain::Domain pmlBnds[3] {std::get<0>(pmlBnd), std::get<1>(pmlBnd), std::get<2>(pmlBnd)}; + for(int i = 0; i < 3; ++i) { + _vC[i] = new gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 >("vC_" + orientation() + "_" + IlearnToCount[i], pmlBnds[i], gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); + } - _vSymC[1][0] = new gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 >("vSymC_second_" + orientation() + "_first", std::get<0>(pmlBnd), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); - _vSymC[1][1] = new gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 >("vSymC_second_" + orientation() + "_second", std::get<1>(pmlBnd), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); - _vSymC[1][2] = new gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 >("vSymC_second_" + orientation() + "_third", std::get<2>(pmlBnd), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); + for(int i = 0; i < 3; ++i) { + for(int j = 0; j < 2; ++j) { + _vSymC[j][i] = new gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 >("vSymC_" + IlearnToCount[j] + "_" + orientation() + "_" + IlearnToCount[i], pmlBnds[i], gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); + } + } - _vSymSymC[0] = new gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 >("vSymSymC_" + orientation() + "_first", std::get<0>(pmlBnd), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); - _vSymSymC[1] = new gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 >("vSymSymC_" + orientation() + "_second", std::get<1>(pmlBnd), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); - _vSymSymC[2] = new gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 >("vSymSymC_" + orientation() + "_third", std::get<2>(pmlBnd), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); + for(int i = 0; i < 3; ++i) { + _vSymSymC[i] = new gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 >("vSymSymC_" + orientation() + "_" + IlearnToCount[i], pmlBnds[i], gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); + } gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 > *uPml[3]; - uPml[0] = static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_bnd[0])->getUPml(); - uPml[1] = static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_bnd[1])->getUPml(); - uPml[2] = static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_bnd[2])->getUPml(); + for(int i = 0; i < 3; ++i) { + uPml[i] = static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_bnd[i])->getUPml(); + } gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 > *uSymPml[2][3]; - uSymPml[0][0] = static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_bnd[0])->getUSymPml(0); - uSymPml[0][1] = static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_bnd[1])->getUSymPml(0); - uSymPml[0][2] = static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_bnd[2])->getUSymPml(0); - - uSymPml[1][0] = static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_bnd[0])->getUSymPml(1); - uSymPml[1][1] = static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_bnd[1])->getUSymPml(1); - uSymPml[1][2] = static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_bnd[2])->getUSymPml(1); + for(int i = 0; i < 3; ++i) { + for(int j = 0; j < 2; ++j) { + uSymPml[j][i] = static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_bnd[i])->getUSymPml(j); + } + } gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 > *uSymSymPml[3]; - uSymSymPml[0] = static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_bnd[0])->getUSymSymPml(); - uSymSymPml[1] = static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_bnd[1])->getUSymSymPml(); - uSymSymPml[2] = static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_bnd[2])->getUSymSymPml(); - + for(int i = 0; i < 3; ++i) { + uSymSymPml[i] = static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_bnd[i])->getUSymSymPml(); + } + std::vector< PmlContinuousSymmetric * > boundaries(6); for(unsigned int i = 0; i < 3; ++i) { boundaries[2 * i] = static_cast< PmlContinuousSymmetric_PmlContinuousSymmetric * >(_bnd[i])->firstBoundary(); @@ -1404,8 +1401,8 @@ namespace D3 { gmshfem::function::TensorFunction< std::complex< double >, 4 > D[3]; gmshfem::function::ScalarFunction< std::complex< double > > E[3]; gmshfem::function::TensorFunction< std::complex< double >, 4 > CMod = boundaries[0]->pmlCoefficients(D[0], E[0], corner, C, kappaP, kappaS); - boundaries[2]->pmlCoefficients(D[1], E[1], corner, C, kappaP, kappaS); - boundaries[4]->pmlCoefficients(D[2], E[2], corner, C, kappaP, kappaS); + boundaries[2]->pmlCoefficients(D[1], E[1], corner, C, fake, fake); + boundaries[4]->pmlCoefficients(D[2], E[2], corner, C, fake, fake); gmshfem::function::ScalarFunction< std::complex< double > > Et = E[0] * E[1] * E[2]; gmshfem::function::VectorFunction< std::complex< double > > n = gmshfem::function::normal< std::complex< double > >(); @@ -1414,9 +1411,7 @@ namespace D3 { gmshfem::function::TensorFunction< std::complex< double > > It = I - In; const std::complex< double > im(0., 1.); - - gmshfem::domain::Domain pmlBnds[3] {std::get<0>(pmlBnd), std::get<1>(pmlBnd), std::get<2>(pmlBnd)}; - + for(unsigned int i = 0; i < 3; ++i) { formulation.integral(- dof(*_vC[i]), tf(*uPml[i]), pmlBnds[i], gauss); @@ -1852,23 +1847,23 @@ namespace D3 { formulation.integral(dof(*_vSymSymCorner[i]), tf(*_vSymSymEdge[i][j]), corner, gauss); } } - - - + + + // third order symmetric fields (symSymSym), u)xyz formulation.integral(- Dt * hxyz * grad(dof(*_uSymSymSymPmlCorner)), hxyz * grad(tf(*_uSymSymSymPmlCorner)), pmlCorner, gauss); formulation.integral(Et * dof(*_uSymSymSymPmlCorner), tf(*_uSymSymSymPmlCorner), pmlCorner, gauss); - + formulation.integral((im/kappaP) * In * dof(*_uSymSymSymPmlCorner), tf(*_uSymSymSymPmlCorner), pmlCornerInf, gauss); formulation.integral((im/kappaS) * It * dof(*_uSymSymSymPmlCorner), tf(*_uSymSymSymPmlCorner), pmlCornerInf, gauss); - + for(int i = 0; i < 3; ++i) { formulation.integral(dof(*_vSymSymSymC[i]), tf(*_uSymSymSymPmlCorner), pmlBnds[i], gauss); formulation.integral(dof(*_uSymSymSymPmlCorner), tf(*_vSymSymSymC[i]), pmlBnds[i], gauss); formulation.integral(- dof(*uSymSymPml[i]), tf(*_vSymSymSymC[i]), pmlBnds[i], gauss); } - + // corner equation for(int i = 0; i < 3; ++i) { _vSymSymSymEdge[i] = new gmshfem::field::CompoundField< std::complex< double >, gmshfem::field::Form::Form0, 3 >("vSymSymSymCornerEdge_" + orientation() + "_" + IlearnToCount[i], cornerEdges[i], gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder); diff --git a/examples/navier/crossPoints/ddm3D.cpp b/examples/navier/crossPoints/ddm3D.cpp index 915e02c0..e4e7486d 100755 --- a/examples/navier/crossPoints/ddm3D.cpp +++ b/examples/navier/crossPoints/ddm3D.cpp @@ -794,19 +794,19 @@ namespace D3 { std::string bnd = boundary; if(boundary == "pml") { if(pmlMethod == "continuous") { - bnd += "Continuous_" + std::to_string(pmlSize) + "_" + pmlType + "_int"; + bnd += "Continuous_" + std::to_string(pmlSize) + "_" + pmlType; } else if(pmlMethod == "continuousSymmetric") { - bnd += "ContinuousSymmetric_" + std::to_string(pmlSize) + "_" + pmlType + "_int"; + bnd += "ContinuousSymmetric_" + std::to_string(pmlSize) + "_" + pmlType; } } std::string bndExt = boundaryExt; if(boundaryExt == "pml") { if(pmlMethodExt == "continuous") { - bndExt += "Continuous_" + std::to_string(pmlSizeExt) + "_" + pmlTypeExt + "_ext"; + bndExt += "Continuous_" + std::to_string(pmlSizeExt) + "_" + pmlTypeExt; } else if(pmlMethodExt == "continuousSymmetric") { - bndExt += "ContinuousSymmetric_" + std::to_string(pmlSizeExt) + "_" + pmlTypeExt + "_ext"; + bndExt += "ContinuousSymmetric_" + std::to_string(pmlSizeExt) + "_" + pmlTypeExt; } } if(gmshddm::mpi::getMPIRank() == 0) { diff --git a/examples/navier/crossPoints/mesh.cpp b/examples/navier/crossPoints/mesh.cpp index cd5522b1..996c2be8 100755 --- a/examples/navier/crossPoints/mesh.cpp +++ b/examples/navier/crossPoints/mesh.cpp @@ -448,7 +448,6 @@ namespace D3 { cubeSurface[3] = newSurface(cubeLines[3], cubeLines[4], cubeLines[11], cubeLines[7]); cubeSurface[4] = newSurface(-cubeLines[0], cubeLines[1], cubeLines[2], cubeLines[3]); cubeSurface[5] = newSurface(cubeLines[8], cubeLines[9], cubeLines[10], cubeLines[11]); - gmsh::model::geo::synchronize(); // gmshfem::msg::warning << "Surfaces : " << cubeSurface[0] << " " << cubeSurface[1] << " " << cubeSurface[2] << " " << cubeSurface[3] << " " << cubeSurface[4] << " " << cubeSurface[5] << gmshfem::msg::endl; int cubeSurfaceLoop = gmsh::model::geo::addSurfaceLoop(cubeSurface); @@ -666,7 +665,7 @@ namespace D3 { { {1, 0, 0, 1, 1, 1}, {0, 0, 0, 0, 1, 1}, - {0, 0, 0, 1, 0, 1} + {0, 1, 0, 1, 1, 1} }, { {0, 1, 0, 1, 1, 1}, @@ -676,7 +675,7 @@ namespace D3 { { {0, 0, 0, 0, 1, 1}, {1, 0, 0, 1, 1, 1}, - {0, 1, 0, 1, 1, 1} + {0, 0, 0, 1, 0, 1} } }; -- GitLab