"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Functions/F_DiffGeom.cpp" between
getdp-3.4.0-source.tgz and getdp-3.5.0-source.tgz

About: GetDP is a general finite element solver using mixed elements to discretize de Rham-type complexes in one, two and three dimensions.

[ To the main getdp source changes report ]

F_DiffGeom.cpp  (getdp-3.4.0-source.tgz):F_DiffGeom.cpp  (getdp-3.5.0-source.tgz)
// GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege // GetDP - Copyright (C) 1997-2022 P. Dular and C. Geuzaine, University of Liege
// //
// See the LICENSE.txt file for license information. Please report all // See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/getdp/getdp/issues. // issues on https://gitlab.onelab.info/getdp/getdp/issues.
// //
// Contributed by Matti Pellikka <matti.pellikka@tut.fi>. // Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
#include <math.h> #include <math.h>
#include "GetDPConfig.h" #include "GetDPConfig.h"
#include "ProData.h" #include "ProData.h"
#include "F.h" #include "F.h"
skipping to change at line 34 skipping to change at line 34
// Arguments: // Arguments:
// k-form coefficient vector // k-form coefficient vector
// metric tensor coefficient matrix // metric tensor coefficient matrix
// //
// Parameters: // Parameters:
// degree k of the form // degree k of the form
// //
// Output: // Output:
// (3-k)-form coefficient vector // (3-k)-form coefficient vector
// //
void F_Hodge(F_ARG) void F_Hodge(F_ARG)
{ {
int k; int k;
struct Value detS; struct Value detS;
struct Value *S; struct Value *S;
k = Fct->Para[0]; k = Fct->Para[0];
S = A+1; S = A + 1;
if( (A->Type != SCALAR && A->Type != VECTOR) || if((A->Type != SCALAR && A->Type != VECTOR) ||
(S->Type != TENSOR_DIAG && (S->Type != TENSOR_DIAG && S->Type != TENSOR_SYM && S->Type != TENSOR))
S->Type != TENSOR_SYM &&
S->Type != TENSOR))
Message::Error("Wrong type of arguments for function 'Hodge'"); Message::Error("Wrong type of arguments for function 'Hodge'");
Cal_DetValue(S, &detS); Cal_DetValue(S, &detS);
detS.Val[0] = sqrt(fabs(detS.Val[0])); detS.Val[0] = sqrt(fabs(detS.Val[0]));
switch(k) { switch(k) {
case 0: case 0: Cal_ProductValue(&detS, A, V); break;
Cal_ProductValue(&detS, A, V);
break;
case 1: case 1:
Cal_InvertValue(S, S); Cal_InvertValue(S, S);
Cal_ProductValue(S, A, V); Cal_ProductValue(S, A, V);
Cal_ProductValue(&detS, V, V); Cal_ProductValue(&detS, V, V);
break; break;
case 2: case 2:
Cal_InvertValue(&detS, &detS); Cal_InvertValue(&detS, &detS);
Cal_ProductValue(S, A, V); Cal_ProductValue(S, A, V);
Cal_ProductValue(&detS, V, V); Cal_ProductValue(&detS, V, V);
break; break;
case 3: case 3:
Cal_InvertValue(&detS, &detS); Cal_InvertValue(&detS, &detS);
Cal_ProductValue(&detS, A, V); Cal_ProductValue(&detS, A, V);
break; break;
default: default: Message::Error("Invalid parameter for function 'Hodge'"); break;
Message::Error("Invalid parameter for function 'Hodge'");
break;
} }
} }
// Inner product of k-forms: // Inner product of k-forms:
// Arguments: // Arguments:
// k-form coefficient vector // k-form coefficient vector
// k-form coefficient vector // k-form coefficient vector
// metric tensor coefficient matrix // metric tensor coefficient matrix
// //
// Parameters: // Parameters:
// degree k of the forms // degree k of the forms
// //
// Output: // Output:
// scalar // scalar
// //
void F_InnerProduct(F_ARG) void F_InnerProduct(F_ARG)
{ {
int k; int k;
struct Value detS; struct Value detS;
struct Value *S; struct Value *S;
struct Value *V1; struct Value *V1;
struct Value *V2; struct Value *V2;
k = Fct->Para[0]; k = Fct->Para[0];
V1 = A; V1 = A;
V2 = A+1; V2 = A + 1;
S = A+2; S = A + 2;
if( (V1->Type != SCALAR && V1->Type != VECTOR) || if((V1->Type != SCALAR && V1->Type != VECTOR) ||
(V2->Type != SCALAR && V2->Type != VECTOR) || (V2->Type != SCALAR && V2->Type != VECTOR) || (V2->Type != V1->Type) ||
(V2->Type != V1->Type) || (S->Type != TENSOR_DIAG && S->Type != TENSOR_SYM && S->Type != TENSOR))
(S->Type != TENSOR_DIAG &&
S->Type != TENSOR_SYM &&
S->Type != TENSOR))
Message::Error("Wrong type of arguments for function 'InnerProduct'"); Message::Error("Wrong type of arguments for function 'InnerProduct'");
switch(k) { switch(k) {
case 0: case 0: Cal_CopyValue(V2, V); break;
Cal_CopyValue(V2, V);
break;
case 1: case 1:
Cal_InvertValue(S, S); Cal_InvertValue(S, S);
Cal_ProductValue(S, V2, V); Cal_ProductValue(S, V2, V);
break; break;
case 2: case 2:
Cal_InvertValue(&detS, &detS); Cal_InvertValue(&detS, &detS);
Cal_ProductValue(S, V2, V); Cal_ProductValue(S, V2, V);
Cal_ProductValue(&detS, V, V); Cal_ProductValue(&detS, V, V);
break; break;
case 3: case 3:
skipping to change at line 145 skipping to change at line 134
// Arguments: // Arguments:
// k-form coefficient vector // k-form coefficient vector
// metric tensor coefficient matrix // metric tensor coefficient matrix
// //
// Parameters: // Parameters:
// degree k of the form // degree k of the form
// //
// Output: // Output:
// k-vector coefficient vector // k-vector coefficient vector
// //
void F_Sharp(F_ARG) void F_Sharp(F_ARG)
{ {
if( (A->Type != SCALAR && A->Type != VECTOR) || if((A->Type != SCALAR && A->Type != VECTOR) ||
((A+1)->Type != TENSOR_DIAG && ((A + 1)->Type != TENSOR_DIAG && (A + 1)->Type != TENSOR_SYM &&
(A+1)->Type != TENSOR_SYM && (A + 1)->Type != TENSOR))
(A+1)->Type != TENSOR))
Message::Error("Wrong type of arguments for function 'Sharp'"); Message::Error("Wrong type of arguments for function 'Sharp'");
if( Fct->Para[0] > 3 || Fct->Para[0] < 0 ) if(Fct->Para[0] > 3 || Fct->Para[0] < 0)
Message::Error("Invalid parameter for function 'Sharp'"); Message::Error("Invalid parameter for function 'Sharp'");
Cal_InvertValue(A+1, A+1); Cal_InvertValue(A + 1, A + 1);
F_Flat(Fct, A, V); F_Flat(Fct, A, V);
} }
// Flat operator of a k-vector: // Flat operator of a k-vector:
// Arguments: // Arguments:
// k-vector coefficient vector // k-vector coefficient vector
// metric tensor coefficient matrix // metric tensor coefficient matrix
// //
// Parameters: // Parameters:
// degree k of the vector // degree k of the vector
// //
// Output: // Output:
// k-form coefficient vector // k-form coefficient vector
// //
void F_Flat(F_ARG) void F_Flat(F_ARG)
{ {
int k; int k;
struct Value detS; struct Value detS;
struct Value *S; struct Value *S;
k = Fct->Para[0]; k = Fct->Para[0];
S = A+1; S = A + 1;
if( (A->Type != SCALAR && A->Type != VECTOR) || if((A->Type != SCALAR && A->Type != VECTOR) ||
(S->Type != TENSOR_DIAG && (S->Type != TENSOR_DIAG && S->Type != TENSOR_SYM && S->Type != TENSOR))
S->Type != TENSOR_SYM &&
S->Type != TENSOR))
Message::Error("Wrong type of arguments for function 'Flat'"); Message::Error("Wrong type of arguments for function 'Flat'");
Cal_DetValue(S, &detS); Cal_DetValue(S, &detS);
detS.Val[0] = sqrt(fabs(detS.Val[0])); detS.Val[0] = sqrt(fabs(detS.Val[0]));
switch(k) { switch(k) {
case 0: case 0: Cal_CopyValue(A, V); break;
Cal_CopyValue(A, V); case 1: Cal_ProductValue(S, A, V); break;
break;
case 1:
Cal_ProductValue(S, A, V);
break;
case 2: case 2:
Cal_InvertValue(S, S); Cal_InvertValue(S, S);
Cal_ProductValue(S, A, V); Cal_ProductValue(S, A, V);
Cal_ProductValue(&detS, V, V); Cal_ProductValue(&detS, V, V);
break; break;
case 3: case 3: Cal_ProductValue(&detS, A, V); break;
Cal_ProductValue(&detS, A, V); default: Message::Error("Invalid parameter for function 'Flat'"); break;
break;
default:
Message::Error("Invalid parameter for function 'Flat'");
break;
} }
} }
// Wedge product of k-forms or k-vectors: // Wedge product of k-forms or k-vectors:
// Arguments: // Arguments:
// k1-form or k1-vector coefficient vector // k1-form or k1-vector coefficient vector
// k2-form or k2-vector coefficient vector // k2-form or k2-vector coefficient vector
// //
// Parameters: // Parameters:
// degree k1 of the first argument // degree k1 of the first argument
// degree k2 of the second argument // degree k2 of the second argument
// //
// Output: // Output:
// (k1+k2)-form or (k1+k2)-vector coefficient vector // (k1+k2)-form or (k1+k2)-vector coefficient vector
// //
void F_WedgeProduct(F_ARG) void F_WedgeProduct(F_ARG)
{ {
int k1,k2; int k1, k2;
struct Value *V1; struct Value *V1;
struct Value *V2; struct Value *V2;
k1 = Fct->Para[0]; k1 = Fct->Para[0];
k2 = Fct->Para[0]; k2 = Fct->Para[0];
V1 = A; V1 = A;
V2 = A+1; V2 = A + 1;
if( (V1->Type != SCALAR && V1->Type != VECTOR) || if((V1->Type != SCALAR && V1->Type != VECTOR) ||
(V2->Type != SCALAR && V2->Type != VECTOR) ) (V2->Type != SCALAR && V2->Type != VECTOR))
Message::Error("Wrong type of arguments for function 'WedgeProduct'"); Message::Error("Wrong type of arguments for function 'WedgeProduct'");
if(k1 < 0 || k1 > 3 || k2 < 0 || k2 > 3) if(k1 < 0 || k1 > 3 || k2 < 0 || k2 > 3)
Message::Error("Invalid parameter for function 'WedgeProduct'"); Message::Error("Invalid parameter for function 'WedgeProduct'");
if( k1 == 0 || k2 == 0 || (k1 == 1 && k2 == 2) || (k1 == 2 && k2 == 1) ) if(k1 == 0 || k2 == 0 || (k1 == 1 && k2 == 2) || (k1 == 2 && k2 == 1))
Cal_ProductValue(V1, V2, V); Cal_ProductValue(V1, V2, V);
else if( k1 == 1 && k2 == 1 ) else if(k1 == 1 && k2 == 1)
Cal_CrossProductValue(V1, V2, V); Cal_CrossProductValue(V1, V2, V);
else if( k1 + k2 > 3 ) else if(k1 + k2 > 3)
Cal_ZeroValue(V); Cal_ZeroValue(V);
else else
Message::Error("Missing implementation in 'WedgeProduct'"); Message::Error("Missing implementation in 'WedgeProduct'");
} }
void F_TensorProduct(F_ARG) void F_TensorProduct(F_ARG)
{ {
Message::Error("'TensorProduct' not implemented"); Message::Error("'TensorProduct' not implemented");
} }
// Interior product of a 1-vector and a k-form: // Interior product of a 1-vector and a k-form:
// Arguments: // Arguments:
// 1-vector coefficient vector // 1-vector coefficient vector
// k-form coefficient vector // k-form coefficient vector
// //
// Parameters: // Parameters:
// degree k of the k-form // degree k of the k-form
// //
// Output: // Output:
// (k-1)-form coefficient vector // (k-1)-form coefficient vector
// //
void F_InteriorProduct(F_ARG) void F_InteriorProduct(F_ARG)
{ {
int k; int k;
struct Value *V1; struct Value *V1;
struct Value *V2; struct Value *V2;
k = Fct->Para[0]; k = Fct->Para[0];
V1 = A; V1 = A;
V2 = A+1; V2 = A + 1;
if( V1->Type != VECTOR || if(V1->Type != VECTOR || (V2->Type != SCALAR && V2->Type != VECTOR))
(V2->Type != SCALAR && V2->Type != VECTOR) )
Message::Error("Wrong type of arguments for function 'InteriorProduct'"); Message::Error("Wrong type of arguments for function 'InteriorProduct'");
switch(k) { switch(k) {
case 1: case 1:
case 3: case 3: Cal_ProductValue(V1, V2, V); break;
Cal_ProductValue(V1, V2, V); case 2: Cal_CrossProductValue(V1, V2, V); break;
break;
case 2:
Cal_CrossProductValue(V1, V2, V);
break;
default: default:
Message::Error("Invalid parameter for function 'InteriorProduct'"); Message::Error("Invalid parameter for function 'InteriorProduct'");
break; break;
} }
} }
// Pullback of a k-form: // Pullback of a k-form:
// Arguments: // Arguments:
// k-form coefficient vector // k-form coefficient vector
// Jacobian matrix of the transition map between charts of a manifold // Jacobian matrix of the transition map between charts of a manifold
// //
// Parameters: // Parameters:
// degree k of the form // degree k of the form
// //
// Output: // Output:
// k-form coefficient vector // k-form coefficient vector
// //
void F_PullBack(F_ARG) void F_PullBack(F_ARG)
{ {
if( (A->Type != SCALAR && A->Type != VECTOR) || if((A->Type != SCALAR && A->Type != VECTOR) ||
((A+1)->Type != TENSOR_DIAG && ((A + 1)->Type != TENSOR_DIAG && (A + 1)->Type != TENSOR_SYM &&
(A+1)->Type != TENSOR_SYM && (A + 1)->Type != TENSOR))
(A+1)->Type != TENSOR))
Message::Error("Wrong type of arguments for function 'PullBack'"); Message::Error("Wrong type of arguments for function 'PullBack'");
if( Fct->Para[0] < 0 || Fct->Para[0] > 3 ) if(Fct->Para[0] < 0 || Fct->Para[0] > 3)
Message::Error("Invalid parameter for function 'PullBack'"); Message::Error("Invalid parameter for function 'PullBack'");
Cal_TransposeValue(A+1, A+1); Cal_TransposeValue(A + 1, A + 1);
F_PushForward(Fct, A, V); F_PushForward(Fct, A, V);
} }
// Pullback of a metric tensor: // Pullback of a metric tensor:
// Arguments: // Arguments:
// metric tensor coefficient matrix // metric tensor coefficient matrix
// Jacobian matrix of the transition map between charts of a manifold // Jacobian matrix of the transition map between charts of a manifold
// //
// Parameters: // Parameters:
// none // none
// //
// Output: // Output:
// metric tensor coefficient matrix // metric tensor coefficient matrix
// //
void F_PullBackMetric(F_ARG) void F_PullBackMetric(F_ARG)
{ {
struct Value *S; struct Value *S;
struct Value *J; struct Value *J;
J = A+1; J = A + 1;
S = A; S = A;
if( (S->Type != TENSOR_DIAG && if((S->Type != TENSOR_DIAG && S->Type != TENSOR_SYM && S->Type != TENSOR) ||
S->Type != TENSOR_SYM && (J->Type != TENSOR_DIAG && J->Type != TENSOR_SYM && J->Type != TENSOR))
S->Type != TENSOR) ||
(J->Type != TENSOR_DIAG &&
J->Type != TENSOR_SYM &&
J->Type != TENSOR))
Message::Error("Wrong type of arguments for function 'PullBackMetric'"); Message::Error("Wrong type of arguments for function 'PullBackMetric'");
Cal_ProductValue(S, J, V); Cal_ProductValue(S, J, V);
Cal_TransposeValue(J, J); Cal_TransposeValue(J, J);
Cal_ProductValue(J, V, V); Cal_ProductValue(J, V, V);
} }
// Inverse pullback of a k-form: // Inverse pullback of a k-form:
// Arguments: // Arguments:
// k-form coefficient vector // k-form coefficient vector
// Jacobian matrix of the transition map between charts of a manifold // Jacobian matrix of the transition map between charts of a manifold
// //
// Parameters: // Parameters:
// degree k of the form // degree k of the form
// //
// Output: // Output:
// k-form coefficient vector // k-form coefficient vector
// //
void F_InvPullBack(F_ARG) void F_InvPullBack(F_ARG)
{ {
if( (A->Type != SCALAR && A->Type != VECTOR) || if((A->Type != SCALAR && A->Type != VECTOR) ||
((A+1)->Type != TENSOR_DIAG && ((A + 1)->Type != TENSOR_DIAG && (A + 1)->Type != TENSOR_SYM &&
(A+1)->Type != TENSOR_SYM && (A + 1)->Type != TENSOR))
(A+1)->Type != TENSOR))
Message::Error("Wrong type of arguments for function 'InvPullBack'"); Message::Error("Wrong type of arguments for function 'InvPullBack'");
if( Fct->Para[0] < 0 || Fct->Para[0] > 3 ) if(Fct->Para[0] < 0 || Fct->Para[0] > 3)
Message::Error("Invalid parameter for function 'InvPullBack'"); Message::Error("Invalid parameter for function 'InvPullBack'");
Cal_InvertValue(A+1, A+1); Cal_InvertValue(A + 1, A + 1);
Cal_TransposeValue(A+1, A+1); Cal_TransposeValue(A + 1, A + 1);
F_PushForward(Fct, A, V); F_PushForward(Fct, A, V);
} }
// Pushforward of a k-vector: // Pushforward of a k-vector:
// Arguments: // Arguments:
// k-vector coefficient vector // k-vector coefficient vector
// Jacobian matrix of the transition map between charts of a manifold // Jacobian matrix of the transition map between charts of a manifold
// //
// Parameters: // Parameters:
// degree k of the k-vector // degree k of the k-vector
// //
// Output: // Output:
// k-vector coefficient vector // k-vector coefficient vector
// //
void F_PushForward(F_ARG) void F_PushForward(F_ARG)
{ {
int k; int k;
struct Value *J; struct Value *J;
struct Value detJ; struct Value detJ;
k = Fct->Para[0]; k = Fct->Para[0];
J = A+1; J = A + 1;
if( (A->Type != SCALAR && A->Type != VECTOR) || if((A->Type != SCALAR && A->Type != VECTOR) ||
(J->Type != TENSOR_DIAG && (J->Type != TENSOR_DIAG && J->Type != TENSOR_SYM && J->Type != TENSOR))
J->Type != TENSOR_SYM &&
J->Type != TENSOR))
Message::Error("Wrong type of arguments for function 'PushForward'"); Message::Error("Wrong type of arguments for function 'PushForward'");
switch(k) { switch(k) {
case 0: case 0: Cal_CopyValue(A, V); break;
Cal_CopyValue(A, V); case 1: Cal_ProductValue(J, A, V); break;
break;
case 1:
Cal_ProductValue(J, A, V);
break;
case 2: case 2:
Cal_InvertValue(J, J); Cal_InvertValue(J, J);
Cal_TransposeValue(J, J); Cal_TransposeValue(J, J);
Cal_DetValue(J, &detJ); Cal_DetValue(J, &detJ);
Cal_ProductValue(J, A, V); Cal_ProductValue(J, A, V);
Cal_ProductValue(&detJ, V, V); Cal_ProductValue(&detJ, V, V);
break; break;
case 3: case 3:
Cal_DetValue(J, &detJ); Cal_DetValue(J, &detJ);
Cal_ProductValue(&detJ, A, V); Cal_ProductValue(&detJ, A, V);
skipping to change at line 441 skipping to change at line 402
// Arguments: // Arguments:
// k-vector coefficient vector // k-vector coefficient vector
// Jacobian matrix of the transition map between charts of a manifold // Jacobian matrix of the transition map between charts of a manifold
// //
// Parameters: // Parameters:
// degree k of the k-vector // degree k of the k-vector
// //
// Output: // Output:
// k-vector coefficient vector // k-vector coefficient vector
// //
void F_InvPushForward(F_ARG) void F_InvPushForward(F_ARG)
{ {
if( (A->Type != SCALAR && A->Type != VECTOR) || if((A->Type != SCALAR && A->Type != VECTOR) ||
((A+1)->Type != TENSOR_DIAG && ((A + 1)->Type != TENSOR_DIAG && (A + 1)->Type != TENSOR_SYM &&
(A+1)->Type != TENSOR_SYM && (A + 1)->Type != TENSOR))
(A+1)->Type != TENSOR))
Message::Error("Wrong type of arguments for function 'InvPushForward'"); Message::Error("Wrong type of arguments for function 'InvPushForward'");
if( Fct->Para[0] < 0 || Fct->Para[0] > 3 ) if(Fct->Para[0] < 0 || Fct->Para[0] > 3)
Message::Error("Invalid parameter for function 'InvPushForward'"); Message::Error("Invalid parameter for function 'InvPushForward'");
Cal_InvertValue(A+1, A+1); Cal_InvertValue(A + 1, A + 1);
F_PushForward(Fct, A, V); F_PushForward(Fct, A, V);
} }
 End of changes. 54 change blocks. 
115 lines changed or deleted 75 lines changed or added
To the C++ Programs section of the getdp source changes report or the top of this page
Mainly generated by pkgdiff using rfcdiff
Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)