diff --git a/examples/helmholtz/crossPoints/Subproblem3D.cpp b/examples/helmholtz/crossPoints/Subproblem3D.cpp index 0b1f042ed2942ffee23e0ab5d4e62eedb6644bca..0517c841a09b31287198e9b34166c2fd45329669 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 c1c27170e4f31f3c94593a82cc1cdc25ef70c803..5b3fd281d198030853ccbd2e99d6da0702f5beb1 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 d82e5ecb99dc3d46c5c44e67d39e9335292d4b48..60eb8bf70ea12250bbaad434141bc2157014c6fc 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 37197a26f8f731b120da3ec7cc682c73536cdb7e..746c859e612282b1c87c956dc4a6e4a6a0b30b15 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 966d1cbd23ab39edf4fa3f6d8d7062391acfc20f..11d03d4c5c0c640f7ea1708aa366be51f6b91387 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 5b3500ef5ad2343fda3d33f6bcb5480dca11183e..fa1f03bff88265bc60cad878172e9881d596a727 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 915e02c0feab3020726acd870e0ee8d6f350b1ca..e4e7486d83081decf7181cf9166b69b6b9b5ede3 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 cd5522b1c8057f4d9b7a622e0fba124d4643cd4f..996c2be82226c62ffbc5f23cbdb7f21881190e6c 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} } };