Skip to content

Commit

Permalink
Merge pull request #98 from UG4/97-limex-does-not-support-norms-on-lo…
Browse files Browse the repository at this point in the history 
…wer-dimensional-entitites-eg-fractures

Using generalized inverses for lower dimensional entities
  • Loading branch information
@anaegel
anaegel authored Dec 3, 2024
2 parents 9851ffa + 5fc64d7 commit a409524
Show file tree
Hide file tree
Showing 3 changed files with 133 additions and 13 deletions.
21 changes: 21 additions & 0 deletions ugbase/common/math/math_vector_matrix/math_matrix_functions.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -401,6 +401,27 @@ inline typename MathMatrix<3,3,T>::value_type
LeftInverse(MathMatrix<3,3>& mOut, const MathMatrix<3,3>& m);
/// \}

////////////////////////////////////////////////////////////////////////////////
// Generalized-Inverse of Matrix
////////////////////////////////////////////////////////////////////////////////
/**
* Computes the pseudo-Inverse of a MxN Matrix (N!=M) and returns the square root of the
* gram determinate or the inverse of a matrix (M=N) and returns the determinante.
*
* For M<N, Right-Inverse of Matrix
* For M>N, Left-Inverse of Matrix
* For M=N, Inverse of Matrix
*
* @param mOut (pseundo)-Inverse of Matrix
* @param m Matrix
* @return Square root of gram determinate of Matrix or determinate
*/
/// \{
template<size_t N, size_t M, typename T>
inline typename MathMatrix<N,M,T>::value_type
GeneralizedInverse(MathMatrix<N,M,T>& mOut, const MathMatrix<M,N,T>& m);
/// \}

////////////////////////////////////////////////////////////////////////////////
// Trace of Matrix
////////////////////////////////////////////////////////////////////////////////
Expand Down
20 changes: 18 additions & 2 deletions ugbase/common/math/math_vector_matrix/math_matrix_functions_common_impl.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -555,12 +555,11 @@ SqrtGramDeterminant(const MathMatrix<3,3,T>& m)
////////////////////////////////////////////////////////////////////////////////
// Inverse of Matrix
////////////////////////////////////////////////////////////////////////////////

template <size_t N, size_t M, typename T>
inline typename MathMatrix<N,M,T>::value_type
Inverse(MathMatrix<N,M,T>& mOut, const MathMatrix<M,N,T>& m)
{
UG_THROW("Inverse for matrix of size "<<N<<"x"<<M<<" not implemented.");
UG_THROW("Inverse for matrix of size "<<M<<"x"<<N<<" not implemented. You could use GeneralizedInverse for pseudo-Inverse.");
}

template <typename T>
Expand Down Expand Up @@ -744,6 +743,23 @@ template <typename T>
inline typename MathMatrix<3,3,T>::value_type
LeftInverse(MathMatrix<3,3>& mOut, const MathMatrix<3,3>& m){return fabs(Inverse(mOut, m));}

////////////////////////////////////////////////////////////////////////////////
// Generalized-Inverse of Matrix
////////////////////////////////////////////////////////////////////////////////
template<size_t N, size_t M, typename T>
inline typename MathMatrix<N,M,T>::value_type
GeneralizedInverse(MathMatrix<N,M,T>& mOut, const MathMatrix<M,N,T>& m)
{
if(M<N){//UG_LOG("Right Inverse for matrix of size "<<M<<"x"<<N<<".");
return RightInverse(mOut,m);
}

if(M>N){//UG_LOG("Left Inverse for matrix of size "<<M<<"x"<<N<<".");
return LeftInverse(mOut,m);
}
return Inverse(mOut,m);
}

////////////////////////////////////////////////////////////////////////////////
// Trace of Matrix
////////////////////////////////////////////////////////////////////////////////
Expand Down
105 changes: 94 additions & 11 deletions ugbase/lib_disc/function_spaces/integrate.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1294,6 +1294,9 @@ class H1ErrorIntegrand
// get reference object id (i.e. Triangle, Quadrilateral, Tetrahedron, ...)
const ReferenceObjectID roid = pElem->reference_object_id();

DimReferenceMapping<elemDim, worldDim>& map
= ReferenceMappingProvider::get<elemDim, worldDim>(roid, vCornerCoords);

try{
// get trial space
const LocalShapeFunctionSet<elemDim>& rTrialSpace =
Expand Down Expand Up @@ -1344,7 +1347,16 @@ class H1ErrorIntegrand
// compute global gradient
MathVector<worldDim> approxGradIP;
MathMatrix<worldDim, elemDim> JTInv;
Inverse(JTInv, vJT[ip]);
map.jacobian_transposed_inverse(JTInv, vLocIP[ip]); //Inverse(JTInv, vJT[ip]);

#ifdef UG_DEBUG
MathMatrix<worldDim, elemDim> JTInv1;
GeneralizedInverse(JTInv1, vJT[ip]); //old algorithms with Inverse, but Inverse is only defined for NxN Matrix
MathMatrix<worldDim, elemDim> diffJTInv;
MatSubtract(diffJTInv, JTInv, JTInv1);
UG_ASSERT(MatMaxNorm(diffJTInv)<SMALL, "The inverse matrix is not identity to the map jocabian transposed inverse.");
#endif //UG_DEBUG

MatVecMult(approxGradIP, JTInv, locTmp);

// get squared of difference
Expand Down Expand Up @@ -1987,6 +1999,9 @@ class H1SemiIntegrand
const ReferenceObjectID roid = pElem->reference_object_id();
const TGridFunction &gridFct= m_scalarData.grid_function();

DimReferenceMapping<elemDim, worldDim>& map
= ReferenceMappingProvider::get<elemDim, worldDim>(roid, vCornerCoords);

typedef typename weight_type::data_type ipdata_type;

std::vector<ipdata_type> elemWeights(numIP, MathMatrix<worldDim, worldDim>());
Expand Down Expand Up @@ -2061,7 +2076,16 @@ class H1SemiIntegrand
// compute gradient
MathVector<worldDim> approxGradIP;
MathMatrix<worldDim, elemDim> JTInv;
Inverse(JTInv, vJT[ip]);
map.jacobian_transposed_inverse(JTInv, vLocIP[ip]);//Inverse(JTInv, vJT[ip]);

#ifdef UG_DEBUG
MathMatrix<worldDim, elemDim> JTInv1;
GeneralizedInverse(JTInv1, vJT[ip]); //old algorithms with Inverse, but Inverse is only defined for NxN Matrix
MathMatrix<worldDim, elemDim> diffJTInv;
MatSubtract(diffJTInv, JTInv, JTInv1);
UG_ASSERT(MatMaxNorm(diffJTInv)<SMALL, "The inverse matrix is not identity to the map jocabian transposed inverse.");
#endif //UG_DEBUG

MatVecMult(approxGradIP, JTInv, tmpVec);

// get norm squared
Expand Down Expand Up @@ -2237,6 +2261,8 @@ class H1SemiDistIntegrand : public StdIntegrand<number, TGridFunction::dim, H1Se
// get reference Mapping
DimReferenceMapping<elemDim, worldDim>& map
= ReferenceMappingProvider::get<elemDim, worldDim>(coarseROID, vCornerCoarse);
DimReferenceMapping<elemDim, worldDim>& mapF
= ReferenceMappingProvider::get<elemDim, worldDim>(fineROID, vCornerCoords);

std::vector<MathVector<elemDim> > vCoarseLocIP;
vCoarseLocIP.resize(numIP);
Expand Down Expand Up @@ -2318,7 +2344,16 @@ class H1SemiDistIntegrand : public StdIntegrand<number, TGridFunction::dim, H1Se
// compute global gradient
MathVector<worldDim> fineGradIP;
MathMatrix<worldDim, elemDim> fineJTInv;
Inverse(fineJTInv, vJT[ip]);
mapF.jacobian_transposed_inverse(fineJTInv, vFineLocIP[ip]);//Inverse(fineJTInv, vJT[ip]);

#ifdef UG_DEBUG
MathMatrix<worldDim, elemDim> fineJTInv1;
GeneralizedInverse(fineJTInv1, vJT[ip]); //old algorithms with Inverse, but Inverse is only defined for NxN Matrix
MathMatrix<worldDim, elemDim> diffJTInv;
MatSubtract(diffJTInv, fineJTInv, fineJTInv1);
UG_ASSERT(MatMaxNorm(diffJTInv)<SMALL, "The inverse matrix is not identity to the map jocabian transposed inverse.");
#endif //UG_DEBUG

MatVecMult(fineGradIP, fineJTInv, fineLocTmp);

// compute global gradient
Expand Down Expand Up @@ -2459,6 +2494,9 @@ class H1EnergyIntegrand
const ReferenceObjectID roid = pElem->reference_object_id();
const TGridFunction &gridFct= m_scalarData.grid_function();

DimReferenceMapping<elemDim, worldDim>& map
= ReferenceMappingProvider::get<elemDim, worldDim>(roid, vCornerCoords);

typedef typename weight_type::data_type ipdata_type;

std::vector<ipdata_type> elemWeights(numIP, MathMatrix<worldDim, worldDim>());
Expand Down Expand Up @@ -2533,7 +2571,16 @@ class H1EnergyIntegrand
// compute gradient
MathVector<worldDim> approxGradIP;
MathMatrix<worldDim, elemDim> JTInv;
Inverse(JTInv, vJT[ip]);
map.jacobian_transposed_inverse(JTInv, vLocIP[ip]);//Inverse(JTInv, vJT[ip]);

#ifdef UG_DEBUG
MathMatrix<worldDim, elemDim> JTInv1;
GeneralizedInverse(JTInv1, vJT[ip]); //old algorithms with Inverse, but Inverse is only defined for NxN Matrix
MathMatrix<worldDim, elemDim> diffJTInv;
MatSubtract(diffJTInv, JTInv, JTInv1);
UG_ASSERT(MatMaxNorm(diffJTInv)<SMALL, "The inverse matrix is not identity to the map jocabian transposed inverse.");
#endif //UG_DEBUG

MatVecMult(approxGradIP, JTInv, tmpVec);

// get norm squared
Expand Down Expand Up @@ -2710,6 +2757,8 @@ class H1EnergyDistIntegrand
// get reference Mapping
DimReferenceMapping<elemDim, worldDim>& map
= ReferenceMappingProvider::get<elemDim, worldDim>(coarseROID, vCornerCoarse);
DimReferenceMapping<elemDim, worldDim>& mapF
= ReferenceMappingProvider::get<elemDim, worldDim>(fineROID, vCornerCoords);

std::vector<MathVector<elemDim> > vCoarseLocIP;
vCoarseLocIP.resize(numIP);
Expand Down Expand Up @@ -2791,20 +2840,31 @@ class H1EnergyDistIntegrand
// compute global D*gradient
MathVector<worldDim> fineGradIP;
MathMatrix<worldDim, elemDim> fineJTInv;
Inverse(fineJTInv, vJT[ip]);
mapF.jacobian_transposed_inverse(fineJTInv, vFineLocIP[ip]);//Inverse(fineJTInv, vJT[ip]);

#ifdef UG_DEBUG
MathMatrix<worldDim, elemDim> fineJTInv1;
GeneralizedInverse(fineJTInv1, vJT[ip]); // old algorithms with Inverse(), but Inverse is only defined for NxN Matrix.
MathMatrix<worldDim, elemDim> diffJTInv;
MatSubtract(diffJTInv, fineJTInv, fineJTInv1);
UG_ASSERT(MatMaxNorm(diffJTInv)<SMALL, "The inverse matrix is not identity to the map jocabian transposed inverse.");
#endif //UG_DEBUG

MatVecMult(fineGradIP, fineJTInv, fineLocTmp);
MatVecMult(fineLocTmp, elemWeights[ip], fineGradIP);
MathVector<worldDim> fineWorldLocTmp(0.0);
MatVecMult(fineWorldLocTmp, elemWeights[ip], fineGradIP);

// compute global D*gradient
MathVector<worldDim> coarseGradIP;
MathMatrix<worldDim, elemDim> coarseJTInv;
map.jacobian_transposed_inverse(coarseJTInv, vCoarseLocIP[ip]);
MatVecMult(coarseGradIP, coarseJTInv, coarseLocTmp);
MatVecMult(coarseLocTmp, elemWeights[ip], coarseGradIP);
MathVector<worldDim> coarseWorldLocTmp(0.0);
MatVecMult(coarseWorldLocTmp, elemWeights[ip], coarseGradIP);

// get squared difference
vValue[ip] = VecDistanceSq(fineLocTmp, coarseLocTmp);

vValue[ip] = VecDistanceSq(fineWorldLocTmp, coarseWorldLocTmp);
//vValue[ip] = VecDistanceSq(fineGradIP, coarseGradIP, elemWeights[ip]);
}

}
Expand Down Expand Up @@ -2893,6 +2953,9 @@ class H1NormIntegrand
// get reference object id (i.e. Triangle, Quadrilateral, Tetrahedron, ...)
const ReferenceObjectID roid = pElem->reference_object_id();

DimReferenceMapping<elemDim, worldDim>& map
= ReferenceMappingProvider::get<elemDim, worldDim>(roid, vCornerCoords);

try{
// get trial space
const LocalShapeFunctionSet<elemDim>& rTrialSpace =
Expand Down Expand Up @@ -2936,7 +2999,16 @@ class H1NormIntegrand
// compute global gradient
MathVector<worldDim> approxGradIP;
MathMatrix<worldDim, elemDim> JTInv;
RightInverse(JTInv, vJT[ip]);
map.jacobian_transposed_inverse(JTInv, vLocIP[ip]);//RightInverse(JTInv, vJT[ip]);

#ifdef UG_DEBUG
MathMatrix<worldDim, elemDim> JTInv1;
RightInverse(JTInv1, vJT[ip]); // old algorithms with RightInverse()
MathMatrix<worldDim, elemDim> diffJTInv;
MatSubtract(diffJTInv, JTInv, JTInv1);
UG_ASSERT(MatMaxNorm(diffJTInv)<SMALL, "The inverse matrix is not identity to the map jocabian transposed inverse.");
#endif //UG_DEBUG

MatVecMult(approxGradIP, JTInv, locTmp);

// get squared of difference
Expand Down Expand Up @@ -3076,6 +3148,8 @@ class H1DistIntegrand
// get Reference Mapping
DimReferenceMapping<elemDim, worldDim>& map
= ReferenceMappingProvider::get<elemDim, worldDim>(coarseROID, vCornerCoarse);
DimReferenceMapping<elemDim, worldDim>& mapF
= ReferenceMappingProvider::get<elemDim, worldDim>(fineROID, vCornerCoords);

std::vector<MathVector<elemDim> > vCoarseLocIP;
vCoarseLocIP.resize(numIP);
Expand Down Expand Up @@ -3125,7 +3199,16 @@ class H1DistIntegrand
// compute global gradient
MathVector<worldDim> fineGradIP;
MathMatrix<worldDim, elemDim> fineJTInv;
RightInverse(fineJTInv, vJT[ip]);
mapF.jacobian_transposed_inverse(fineJTInv, vFineLocIP[ip]);//RightInverse(fineJTInv, vJT[ip]);

#ifdef UG_DEBUG
MathMatrix<worldDim, elemDim> fineJTInv1;
RightInverse(fineJTInv1, vJT[ip]); // old algorithms with RightInverse()
MathMatrix<worldDim, elemDim> diffJTInv;
MatSubtract(diffJTInv, fineJTInv, fineJTInv1);
UG_ASSERT(MatMaxNorm(diffJTInv)<SMALL, "The inverse matrix is not identity to the map jocabian transposed inverse.");
#endif //UG_DEBUG

MatVecMult(fineGradIP, fineJTInv, fineLocTmp);

// compute global gradient
Expand Down

0 comments on commit a409524

Please sign in to comment.