diff --git a/dumux/assembly/boxlocalresidual.hh b/dumux/assembly/boxlocalresidual.hh index eaadcedddf6ae9fd64388e42bb5e31f84f77775b..9f105cabc40ed8c3603efc35ede617f0fdb4f428 100644 --- a/dumux/assembly/boxlocalresidual.hh +++ b/dumux/assembly/boxlocalresidual.hh @@ -106,23 +106,21 @@ public: const auto& scv = fvGeometry.scv(scvf.insideScvIdx()); const auto& bcTypes = elemBcTypes[scv.localDofIndex()]; - // Neumann and Robin ("solution dependent Neumann") boundary conditions - if (bcTypes.hasNeumann() && !bcTypes.hasDirichlet()) + // Treat Neumann and Robin ("solution dependent Neumann") boundary conditions. + // For Dirichlet there is no addition to the residual here but they + // are enforced strongly by replacing the residual entry afterwards. + if (bcTypes.hasNeumann()) { auto neumannFluxes = problem.neumann(element, fvGeometry, elemVolVars, scvf); // multiply neumann fluxes with the area and the extrusion factor neumannFluxes *= scvf.area()*elemVolVars[scv].extrusionFactor(); - flux += neumannFluxes; + // only add fluxes to equations for which Neumann is set + for (int eqIdx = 0; eqIdx < NumEqVector::dimension; ++eqIdx) + if (bcTypes.isNeumann(eqIdx)) + flux[eqIdx] += neumannFluxes[eqIdx]; } - - // for Dirichlet there is no addition to the residual here but they - // are enforced strongly by replacing the residual entry afterwards - else if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann()) - return flux; - else - DUNE_THROW(Dune::NotImplemented, "Mixed boundary conditions. Use pure boundary conditions by converting Dirichlet BCs to Robin BCs"); } return flux; diff --git a/test/geomechanics/elastic/problem.hh b/test/geomechanics/elastic/problem.hh index 9762336500325cafcaa8139463a4420475def5b0..b8d2a4053cf0bf2a010e592fa70905928f1ca221 100644 --- a/test/geomechanics/elastic/problem.hh +++ b/test/geomechanics/elastic/problem.hh @@ -107,6 +107,39 @@ public: { BoundaryTypes values; values.setAllDirichlet(); + + if (globalPos[0] < eps_ && globalPos[1] > eps_ && globalPos[1] < 1.0 - eps_) + values.setNeumann(Indices::momentum(/*x-dir*/0)); + if (globalPos[1] > 1.0 - eps_ && globalPos[0] > eps_ && globalPos[0] < 1.0 - eps_) + values.setNeumann(Indices::momentum(/*y-dir*/1)); + + return values; + } + + NumEqVector neumann(const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolvars, + const SubControlVolumeFace& scvf) const + { + GradU gradU = exactGradient(scvf.ipGlobal()); + + GradU epsilon; + for (int i = 0; i < dim; ++i) + for (int j = 0; j < dimWorld; ++j) + epsilon[i][j] = 0.5*(gradU[i][j] + gradU[j][i]); + + const auto& lameParams = this->spatialParams().lameParamsAtPos(scvf.ipGlobal()); + GradU sigma(0.0); + const auto traceEpsilon = trace(epsilon); + for (int i = 0; i < dim; ++i) + { + sigma[i][i] = lameParams.lambda()*traceEpsilon; + for (int j = 0; j < dimWorld; ++j) + sigma[i][j] += 2.0*lameParams.mu()*epsilon[i][j]; + } + + NumEqVector values; + sigma.mv(scvf.unitOuterNormal(), values); return values; } diff --git a/test/references/elasticbox-reference.vtu b/test/references/elasticbox-reference.vtu index 4a1cd64dddd8823444f82b8fd791d9c237103fb3..3199fcc9800a2964f648009ebc536bd02df8675a 100644 --- a/test/references/elasticbox-reference.vtu +++ b/test/references/elasticbox-reference.vtu @@ -4,36 +4,36 @@ <Piece NumberOfCells="100" NumberOfPoints="121"> <PointData Vectors="u"> <DataArray type="Float32" Name="u" NumberOfComponents="3" format="ascii"> - 0 0 0 0 0 0 0 0 0 0.0486269 0.363944 0 - 0 0 0 0.0940417 0.589125 0 0 0 0 0.131474 0.588995 0 - 0 0 0 0.156265 0.363961 0 0 0 0 0.164956 9.03744e-18 0 - 0 0 0 0.156265 -0.363961 0 0 0 0 0.131474 -0.588995 0 - 0 0 0 0.0940417 -0.589125 0 0 0 0 0.0486269 -0.363944 0 - 0 0 0 0 0 0 0 0 0 0.084587 0.589637 0 - 0.15322 0.951768 0 0.204761 0.950788 0 0.236895 0.587348 0 0.247828 3.63266e-17 0 - 0.236895 -0.587348 0 0.204761 -0.950788 0 0.15322 -0.951768 0 0.084587 -0.589637 0 - 0 0 0 0 0 0 0.0876187 0.590116 0 0.151736 0.951624 0 - 0.194708 0.950165 0 0.219141 0.586836 0 0.227034 3.83134e-17 0 0.219141 -0.586836 0 - 0.194708 -0.950165 0 0.151736 -0.951624 0 0.0876187 -0.590116 0 0 0 0 - 0 0 0 0.0560947 0.365428 0 0.0900437 0.58885 0 0.106944 0.587692 0 - 0.113462 0.362888 0 0.114957 7.95142e-18 0 0.113462 -0.362888 0 0.106944 -0.587692 0 - 0.0900437 -0.58885 0 0.0560947 -0.365428 0 0 0 0 0 0 0 - 0.00219115 0.00143062 0 -0.00779061 0.00174145 0 -0.0239508 0.00143597 0 -0.0381341 0.000793643 0 - -0.043704 -6.0621e-17 0 -0.0381341 -0.000793643 0 -0.0239508 -0.00143597 0 -0.00779061 -0.00174145 0 - 0.00219115 -0.00143062 0 0 0 0 0 0 0 -0.0534173 -0.362979 0 - -0.104191 -0.585832 0 -0.147642 -0.585182 0 -0.177287 -0.361494 0 -0.187839 -1.27627e-16 0 - -0.177287 0.361494 0 -0.147642 0.585182 0 -0.104191 0.585832 0 -0.0534173 0.362979 0 - 0 0 0 0 0 0 -0.0895704 -0.588716 0 -0.162557 -0.949787 0 - -0.21726 -0.948574 0 -0.251351 -0.585934 0 -0.262955 -1.50929e-16 0 -0.251351 0.585934 0 - -0.21726 0.948574 0 -0.162557 0.949787 0 -0.0895704 0.588716 0 0 0 0 - 0 0 0 -0.0925845 -0.589427 0 -0.161108 -0.951263 0 -0.207287 -0.950244 0 - -0.233697 -0.587011 0 -0.242265 -1.26245e-16 0 -0.233697 0.587011 0 -0.207287 0.950244 0 - -0.161108 0.951263 0 -0.0925845 0.589427 0 0 0 0 0 0 0 - -0.0608212 -0.364415 0 -0.100277 -0.589403 0 -0.123287 -0.58911 0 -0.13469 -0.363994 0 - -0.138052 -4.5718e-17 0 -0.13469 0.363994 0 -0.123287 0.58911 0 -0.100277 0.589403 0 - -0.0608212 0.364415 0 0 0 0 0 0 0 0 0 0 - 0 0 0 0 0 0 0 0 0 0 0 0 - 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 -0.0167977 0 0 0.0409281 0.363294 0 + 0 0 0 0.0911066 0.588056 0 0 0 0 0.130874 0.588076 0 + 0 0 0 0.156692 0.363309 0 0 0 0 0.165728 -0.000415497 0 + 0 0 0 0.157041 -0.364207 0 0 0 0 0.132098 -0.589139 0 + 0 0 0 0.0944587 -0.589224 0 0 0 0 0.0488348 -0.364026 0 + 0 0 0 0 0 0 -0.00961956 0 0 0.0803135 0.58624 0 + 0.152524 0.948628 0 0.205857 0.948544 0 0.238681 0.585951 0 0.249691 -0.00076681 0 + 0.238514 -0.587702 0 0.205986 -0.950914 0 0.154006 -0.95181 0 0.0849515 -0.589675 0 + 0 0 0 0.00730444 0 0 0.0938215 0.585801 0 0.157096 0.947604 0 + 0.19936 0.947384 0 0.223104 0.585208 0 0.230294 -0.000762067 0 0.221677 -0.587025 0 + 0.196515 -0.950047 0 0.152848 -0.951428 0 0.0881119 -0.59 0 0 0 0 + 0.0225539 0 0 0.0722235 0.362433 0 0.101357 0.585937 0 0.115071 0.585665 0 + 0.119493 0.361814 0 0.119495 -0.000262609 0 0.116801 -0.362566 0 0.109235 -0.587066 0 + 0.0914051 -0.588219 0 0.056676 -0.365044 0 0 0 0 0.0282842 0 0 + 0.0221021 0.00107396 0 0.00586941 0.00118698 0 -0.0143626 0.00100342 0 -0.0311378 0.000792992 0 + -0.038495 0.000591673 0 -0.034339 0.000328459 0 -0.0213882 -5.23562e-05 0 -0.00630706 -0.000480133 0 + 0.00279812 -0.000666758 0 0 0 0 0.0217288 0 0 -0.0378753 -0.360865 0 + -0.0931049 -0.584272 0 -0.139371 -0.584245 0 -0.170833 -0.360546 0 -0.182769 0.00146279 0 + -0.173488 0.36356 0 -0.145081 0.58753 0 -0.102756 0.587909 0 -0.0528684 0.364228 0 + 0 0 0 0.00585659 0 0 -0.0844141 -0.585792 0 -0.157583 -0.94787 0 + -0.212287 -0.947578 0 -0.246536 -0.584814 0 -0.258645 0.00202466 0 -0.247931 0.588974 0 + -0.214971 0.952086 0 -0.161368 0.952873 0 -0.0891822 0.590539 0 0 0 0 + -0.0113062 0 0 -0.0981222 -0.588299 0 -0.162233 -0.951533 0 -0.205688 -0.950996 0 + -0.230703 -0.586782 0 -0.238948 0.00216631 0 -0.230896 0.591086 0 -0.205508 0.955218 0 + -0.160375 0.955605 0 -0.0924721 0.591866 0 0 0 0 -0.0181253 0 0 + -0.0696093 -0.367307 0 -0.103544 -0.593653 0 -0.123354 -0.59295 0 -0.133074 -0.365463 0 + -0.13587 0.00205104 0 -0.132799 0.369326 0 -0.122239 0.596046 0 -0.100158 0.595378 0 + -0.0610389 0.367376 0 0 0 0 0 0 0 0 -0.00525857 0 + 0 -0.00773072 0 0 -0.00697102 0 0 -0.00330957 0 0 0.00196536 0 + 0 0.00689845 0 0 0.00944745 0 0 0.00816588 0 0 0.00313597 0 0 0 0 </DataArray> <DataArray type="Float32" Name="u_exact" NumberOfComponents="3" format="ascii">