Commit a762c728 authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

[decoupled] generalize velocity output to handle also simplices

So far, velocity output for the decoupled models has been restricted to
cubes. It now also works for simplices such as triangles and
tetrahedrons. For each element, it evaluates the Raviart-Thomas-0
interpolant of the fluxes in the barycenter.

Brought to attention by, discussed with and reviewed by Nicolas.



git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@13610 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 82d945d5
......@@ -144,22 +144,27 @@ public:
flux[isIndex] = isIt->geometry().volume()
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(isIndex));
}
Dune::FieldVector<Scalar, dim> refVelocity(0);
//2D simplices
if ( dim==2 && geometry.corners() == 3 ) {
refVelocity[0] = 1.0/3.0*(2.0*flux[1] - flux[2]);
refVelocity[1] = 1.0/3.0*(2.0*flux[0] - flux[2]);
}
//3D tetrahedra
else if ( dim==3 && geometry.corners() == 4 ){
DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
// calculate velocity on reference element as the Raviart-Thomas-0
// interpolant of the fluxes
Dune::FieldVector<Scalar, dim> refVelocity;
// simplices
if (refElement.type().isSimplex()) {
for (int dimIdx = 0; dimIdx < dim; dimIdx++)
{
refVelocity[dimIdx] = -flux[dim - 1 - dimIdx];
for (int faceIdx = 0; faceIdx < dim + 1; faceIdx++)
{
refVelocity[dimIdx] += flux[faceIdx]/(dim + 1);
}
}
}
//2D and 3D cubes
else if ( refElement.type().isCube() ){
// cubes
else if (refElement.type().isCube()){
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (flux[2*i + 1] - flux[2*i]);
}
//3D prism and pyramids
// 3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
......
......@@ -220,25 +220,30 @@ public:
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
}
Dune::FieldVector < Scalar, dim > refVelocity(0);
//2D simplices
if ( dim==2 && geometry.corners() == 3 ) {
refVelocity[0] = 1.0/3.0*(2.0*fluxW[1] - fluxW[2]);
refVelocity[1] = 1.0/3.0*(2.0*fluxW[0] - fluxW[2]);
// calculate velocity on reference element as the Raviart-Thomas-0
// interpolant of the fluxes
Dune::FieldVector<Scalar, dim> refVelocity;
// simplices
if (refElement.type().isSimplex()) {
for (int dimIdx = 0; dimIdx < dim; dimIdx++)
{
refVelocity[dimIdx] = -fluxW[dim - 1 - dimIdx];
for (int faceIdx = 0; faceIdx < dim + 1; faceIdx++)
{
refVelocity[dimIdx] += fluxW[faceIdx]/(dim + 1);
}
}
}
//3D tetrahedra
else if ( dim==3 && geometry.corners() == 4 ){
DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
}
//2D and 3D cubes
else if ( refElement.type().isCube() ){
// cubes
else if (refElement.type().isCube()){
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (fluxW[2*i + 1] - fluxW[2*i]);
}
//3D prism and pyramids
// 3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
const Dune::FieldVector<Scalar, dim> localPos =
refElement.position(0, 0);
......@@ -253,25 +258,29 @@ public:
velocity[globalIdx] = elementVelocity;
refVelocity = 0;
//2D simplices
if ( dim==2 && geometry.corners() == 3 ) {
refVelocity[0] = 1.0/3.0*(2.0*fluxNw[1] - fluxNw[2]);
refVelocity[1] = 1.0/3.0*(2.0*fluxNw[0] - fluxNw[2]);
// calculate velocity on reference element as the Raviart-Thomas-0
// interpolant of the fluxes
// simplices
if (refElement.type().isSimplex()) {
for (int dimIdx = 0; dimIdx < dim; dimIdx++)
{
refVelocity[dimIdx] = -fluxNw[dim - 1 - dimIdx];
for (int faceIdx = 0; faceIdx < dim + 1; faceIdx++)
{
refVelocity[dimIdx] += fluxNw[faceIdx]/(dim + 1);
}
}
}
//3D tetrahedra
else if ( dim==3 && geometry.corners() == 4 ){
DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
}
//2D and 3D cubes
else if ( refElement.type().isCube() ){
// cubes
else if (refElement.type().isCube()){
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (fluxNw[2*i + 1] - fluxNw[2*i]);
}
//3D prism and pyramids
// 3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
// calculate the element velocity by the Piola transformation
elementVelocity = 0;
jacobianT.umtv(refVelocity, elementVelocity);
......
......@@ -302,25 +302,30 @@ public:
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
}
DimVector refVelocity(0.0);
//2D simplices
if ( dim==2 && geometry.corners() == 3 ) {
refVelocity[0] = 1.0/3.0*(2.0*fluxW[1] - fluxW[2]);
refVelocity[1] = 1.0/3.0*(2.0*fluxW[0] - fluxW[2]);
// calculate velocity on reference element as the Raviart-Thomas-0
// interpolant of the fluxes
Dune::FieldVector<Scalar, dim> refVelocity;
// simplices
if (refElement.type().isSimplex()) {
for (int dimIdx = 0; dimIdx < dim; dimIdx++)
{
refVelocity[dimIdx] = -fluxW[dim - 1 - dimIdx];
for (int faceIdx = 0; faceIdx < dim + 1; faceIdx++)
{
refVelocity[dimIdx] += fluxW[faceIdx]/(dim + 1);
}
}
}
//3D tetrahedra
else if ( dim==3 && geometry.corners() == 4 ){
DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
// cubes
else if (refElement.type().isCube()){
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (fluxW[2*i + 1] - fluxW[2*i]);
}
//2D and 3D cubes
else if ( refElement.type().isCube() ){
for (int k = 0; k < dim; k++)
refVelocity[k] = 0.5 * (fluxW[2*k+1] - fluxW[2*k]);
}
//3D prism and pyramids
// 3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
const DimVector& localPos = refElement.position(0, 0);
// get the transposed Jacobian of the element mapping
......@@ -333,24 +338,29 @@ public:
(*velocityWetting)[globalIdx] = elementVelocity;
//2D simplices
if ( dim==2 && geometry.corners() == 3 ) {
refVelocity[0] = 1.0/3.0*(2.0*fluxNw[1] - fluxNw[2]);
refVelocity[1] = 1.0/3.0*(2.0*fluxNw[0] - fluxNw[2]);
// calculate velocity on reference element as the Raviart-Thomas-0
// interpolant of the fluxes
// simplices
if (refElement.type().isSimplex()) {
for (int dimIdx = 0; dimIdx < dim; dimIdx++)
{
refVelocity[dimIdx] = -fluxNw[dim - 1 - dimIdx];
for (int faceIdx = 0; faceIdx < dim + 1; faceIdx++)
{
refVelocity[dimIdx] += fluxNw[faceIdx]/(dim + 1);
}
}
}
//3D tetrahedra
else if ( dim==3 && geometry.corners() == 4 ){
DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
// cubes
else if (refElement.type().isCube()){
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (fluxNw[2*i + 1] - fluxNw[2*i]);
}
//2D and 3D cubes
else if ( refElement.type().isCube() ){
for (int k = 0; k < dim; k++)
refVelocity[k] = 0.5 * (fluxNw[2*k+1] - fluxNw[2*k]);
}
//3D prism and pyramids
// 3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
// calculate the element velocity by the Piola transformation
elementVelocity = 0;
jacobianT.umtv(refVelocity, elementVelocity);
......
......@@ -314,25 +314,30 @@ public:
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
}
Dune::FieldVector < Scalar, dim > refVelocity(0);
//2D simplices
if ( dim==2 && geometry.corners() == 3 ) {
refVelocity[0] = 1.0/3.0*(2.0*fluxW[1] - fluxW[2]);
refVelocity[1] = 1.0/3.0*(2.0*fluxW[0] - fluxW[2]);
// calculate velocity on reference element as the Raviart-Thomas-0
// interpolant of the fluxes
Dune::FieldVector<Scalar, dim> refVelocity;
// simplices
if (refElement.type().isSimplex()) {
for (int dimIdx = 0; dimIdx < dim; dimIdx++)
{
refVelocity[dimIdx] = -fluxW[dim - 1 - dimIdx];
for (int faceIdx = 0; faceIdx < dim + 1; faceIdx++)
{
refVelocity[dimIdx] += fluxW[faceIdx]/(dim + 1);
}
}
}
//3D tetrahedra
else if ( dim==3 && geometry.corners() == 4 ){
DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
// cubes
else if (refElement.type().isCube()){
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (fluxW[2*i + 1] - fluxW[2*i]);
}
//2D and 3D cubes
else if ( refElement.type().isCube() ){
for (int k = 0; k < dim; k++)
refVelocity[k] = 0.5 * (fluxW[2*k+1] - fluxW[2*k]);
}
//3D prism and pyramids
// 3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
const DimVector& localPos = refElement.position(0, 0);
// get the transposed Jacobian of the element mapping
......@@ -345,24 +350,29 @@ public:
(*velocityWetting)[globalIdx] = elementVelocity;
//2D simplices
if ( dim==2 && geometry.corners() == 3 ) {
refVelocity[0] = 1.0/3.0*(2.0*fluxNw[1] - fluxNw[2]);
refVelocity[1] = 1.0/3.0*(2.0*fluxNw[0] - fluxNw[2]);
// calculate velocity on reference element as the Raviart-Thomas-0
// interpolant of the fluxes
// simplices
if (refElement.type().isSimplex()) {
for (int dimIdx = 0; dimIdx < dim; dimIdx++)
{
refVelocity[dimIdx] = -fluxNw[dim - 1 - dimIdx];
for (int faceIdx = 0; faceIdx < dim + 1; faceIdx++)
{
refVelocity[dimIdx] += fluxNw[faceIdx]/(dim + 1);
}
}
}
//3D tetrahedra
else if ( dim==3 && geometry.corners() == 4 ){
DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
// cubes
else if (refElement.type().isCube()){
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (fluxNw[2*i + 1] - fluxNw[2*i]);
}
//2D and 3D cubes
else if ( refElement.type().isCube() ){
for (int k = 0; k < dim; k++)
refVelocity[k] = 0.5 * (fluxNw[2*k+1] - fluxNw[2*k]);
}
//3D prism and pyramids
// 3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
// calculate the element velocity by the Piola transformation
elementVelocity = 0;
jacobianT.umtv(refVelocity, elementVelocity);
......
......@@ -338,22 +338,26 @@ public:
if (useFluxInd_ || usePercentileFlux_)
{
DimVector refVelocity(0);
//2D simplices
if ( dim==2 && geometry.corners() == 3 ) {
refVelocity[0] = 1.0/3.0*(2.0*flux[1] - flux[2]);
refVelocity[1] = 1.0/3.0*(2.0*flux[0] - flux[2]);
}
//3D tetrahedra
else if ( dim==3 && geometry.corners() == 4 ){
DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
// calculate velocity on reference element as the Raviart-Thomas-0
// interpolant of the fluxes
Dune::FieldVector<Scalar, dim> refVelocity;
// simplices
if (refElement.type().isSimplex()) {
for (int dimIdx = 0; dimIdx < dim; dimIdx++)
{
refVelocity[dimIdx] = -flux[dim - 1 - dimIdx];
for (int faceIdx = 0; faceIdx < dim + 1; faceIdx++)
{
refVelocity[dimIdx] += flux[faceIdx]/(dim + 1);
}
}
}
//2D and 3D cubes
else if ( refElement.typ().isCube() ){
// cubes
else if (refElement.type().isCube()){
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (flux[2*i + 1] - flux[2*i]);
}
//3D prism and pyramids
// 3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
......
......@@ -200,15 +200,28 @@ public:
}
}
// calculate velocity on reference element
Dune::FieldVector<Scalar,dim> refVelocity;
if (geometry.corners() == 3) {
refVelocity[0] = 1.0/3.0*(fluxVector[0] + fluxVector[2] - 2.0*fluxVector[1]);
refVelocity[1] = 1.0/3.0*(fluxVector[0] + fluxVector[1] - 2.0*fluxVector[2]);
// calculate velocity on reference element as the Raviart-Thomas-0
// interpolant of the fluxes
Dune::FieldVector<Scalar, dim> refVelocity;
// simplices
if (geometry.type().isSimplex()) {
for (int dimIdx = 0; dimIdx < dim; dimIdx++)
{
refVelocity[dimIdx] = -fluxVector[dim - 1 - dimIdx];
for (int faceIdx = 0; faceIdx < dim + 1; faceIdx++)
{
refVelocity[dimIdx] += fluxVector[faceIdx]/(dim + 1);
}
}
}
// cubes
else if (geometry.type().isCube()){
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (fluxVector[2*i + 1] - fluxVector[2*i]);
}
// 3D prism and pyramids
else {
refVelocity[0] = 0.5*(fluxVector[1] - fluxVector[0]);
refVelocity[1] = 0.5*(fluxVector[3] - fluxVector[2]);
DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
// get the transposed Jacobian of the element mapping
......
......@@ -303,15 +303,28 @@ public:
}
}
// calculate velocity on reference element
Dune::FieldVector<ct,dim> refVelocity;
if (geometry.corners() == 8) {
refVelocity[0] = 0.5*(fluxVector[1] - fluxVector[0]);
refVelocity[1] = 0.5*(fluxVector[3] - fluxVector[2]);
refVelocity[2] = 0.5*(fluxVector[5] - fluxVector[4]);
// calculate velocity on reference element as the Raviart-Thomas-0
// interpolant of the fluxes
Dune::FieldVector<double, dim> refVelocity;
// simplices
if (geometry.type().isSimplex()) {
for (int dimIdx = 0; dimIdx < dim; dimIdx++)
{
refVelocity[dimIdx] = -fluxVector[dim - 1 - dimIdx];
for (int faceIdx = 0; faceIdx < dim + 1; faceIdx++)
{
refVelocity[dimIdx] += fluxVector[faceIdx]/(dim + 1);
}
}
}
// cubes
else if (geometry.type().isCube()){
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (fluxVector[2*i + 1] - fluxVector[2*i]);
}
// 3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented, "grid shape is not supported!");
DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
// get the transposed Jacobian of the element mapping
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment