[cig-commits] r14789 - short/3D/PyLith/trunk/libsrc/materials
willic3 at geodynamics.org
willic3 at geodynamics.org
Thu Apr 23 21:52:48 PDT 2009
Author: willic3
Date: 2009-04-23 21:52:48 -0700 (Thu, 23 Apr 2009)
New Revision: 14789
Modified:
short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh
short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.icc
Log:
More work on PowerLaw3D.
Right now, I think it will be easier and more efficient to write our own
1D root finder. The PETSc method seems too complex with too much
overhead for this.
Modified: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2009-04-23 02:32:01 UTC (rev 14788)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2009-04-24 04:52:48 UTC (rev 14789)
@@ -47,7 +47,7 @@
{ "viscosity_coeff", 1, OTHER_FIELD },
{ "power_law_exponent", 1, OTHER_FIELD },
{ "total_strain", tensorSize, OTHER_FIELD },
- { "viscous_strain", tensorSize, OTHER_FIELD },
+ { "viscous_strain_t", tensorSize, OTHER_FIELD },
{ "stress_t", tensorSize, OTHER_FIELD },
};
/// Indices (order) of properties.
@@ -57,8 +57,8 @@
const int pidViscosityCoeff = pidLambda + 1;
const int pidPowerLawExp = pidViscosityCoeff + 1;
const int pidStrainT = pidPowerLawExp + 1;
- const int pidVisStrain = pidStrainT + tensorSize;
- const int pidStressT = pidVisStrain + tensorSize;
+ const int pidVisStrainT = pidStrainT + tensorSize;
+ const int pidStressT = pidVisStrainT + tensorSize;
/// Values expected in spatial database
const int numDBValues = 5;
@@ -97,7 +97,7 @@
_updatePropertiesFn(&pylith::materials::PowerLaw3D::_updatePropertiesElastic)
{ // constructor
_dimension = 3;
- _visStrain.resize(_PowerLaw3D::tensorSize);
+ _visStrainT.resize(_PowerLaw3D::tensorSize);
_stressT.resize(_PowerLaw3D::tensorSize);
} // constructor
@@ -279,7 +279,7 @@
// ----------------------------------------------------------------------
// Compute viscous strain for current time step.
-//********** May not need this. Check formulation. I probably still need
+//********** May not need this. Check formulation. I still need
// updateState. *************
void
pylith::materials::PowerLaw3D::_computeStateVars(
@@ -326,8 +326,8 @@
devStrainTpdt = totalStrain[iComp] - diag[iComp] * meanStrainTpdt;
devStrainT = properties[_PowerLaw3D::pidStrainT+iComp] -
diag[iComp] * meanStrainT;
- _visStrain[iComp] = expFac *
- properties[_PowerLaw3D::pidVisStrain + iComp] +
+ _visStrainT[iComp] = expFac *
+ properties[_PowerLaw3D::pidVisStrainT + iComp] +
dq * (devStrainTpdt - devStrainT);
} // for
@@ -385,25 +385,24 @@
// Effective stress function that computes effective stress function only
// (no derivative).
double
-pylith::materials::PowerLaw3D::_effStressFunc(double x,
- void *params)
+pylith::materials::PowerLaw3D::_effStressFunc(double effStressTpdt,
+ double *params)
{ // _effStressFunc
- struct effStressParams *p = (struct effStressParams *) params;
- double ae = p->ae;
- double b = p->b;
- double c = p->c;
- double d = p->d;
- double alfa = p->alfa;
- double dt = p->dt;
- double effStressT = p->effStressT;
- double powerLawExp = p->powerLawExp;
- double viscosityCoeff = p->viscosityCoeff;
- double factor1 = 1.0-alfa;
- double effStressTau = factor1 * effStressT + alfa * x;
+ double ae = params[0];
+ double b = params[1];
+ double c = params[2];
+ double d = params[3];
+ double alpha = params[4];
+ double dt = params[5];
+ double effStressT = params[6];
+ double powerLawExp = params[7];
+ double viscosityCoeff = params[8];
+ double factor1 = 1.0-alpha;
+ double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
double gammaTau = 0.5 * ((effStressTau/viscosityCoeff)^
(powerLawExp - 1.0))/viscosityCoeff;
- double a = ae + alfa * dt * gammaTau;
- double y = a * a * x * x - b +
+ double a = ae + alpha * dt * gammaTau;
+ double y = a * a * effStressTpdt * effStressTpdt - b +
c * gammaTau - d * d * gammaTau * gammaTau;
PetscLogFlops(21);
return y;
@@ -413,28 +412,28 @@
// Effective stress function that computes effective stress function
// derivative only (no function value).
double
-pylith::materials::PowerLaw3D::_effStressDFunc(double x,
- void *params)
+pylith::materials::PowerLaw3D::_effStressDFunc(double effStressTpdt,
+ double *params)
{ // _effStressDFunc
- struct effStressParams *p = (struct effStressParams *) params;
- double ae = p->ae;
- double c = p->c;
- double d = p->d;
- double alfa = p->alfa;
- double dt = p->dt;
- double effStressT = p->effStressT;
- double powerLawExp = p->powerLawExp;
- double viscosityCoeff = p->viscosityCoeff;
- double factor1 = 1.0-alfa;
- double effStressTau = factor1 * effStressT + alfa * x;
+ double ae = params[0];
+ double c = params[2];
+ double d = params[3];
+ double alpha = params[4];
+ double dt = params[5];
+ double effStressT = params[6];
+ double powerLawExp = params[7];
+ double viscosityCoeff = params[8];
+ double factor1 = 1.0-alpha;
+ double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
double gammaTau = 0.5 * ((effStressTau/viscosityCoeff)^
(powerLawExp - 1.0))/viscosityCoeff;
- double a = ae + alfa * dt * gammaTau;
- double dGammaTau = 0.5 * alfa * (powerLawExp - 1.0) *
+ double a = ae + alpha * dt * gammaTau;
+ double dGammaTau = 0.5 * alpha * (powerLawExp - 1.0) *
((effStressTau/viscosityCoeff) ^ (powerLawExp - 2.0))/
(viscosityCoeff * viscosityCoeff);
- double dy = 2.0 * a * a * x + dGammaTau *
- (2.0 * a * alfa * dt * x * x + c - 2.0 * d * d * gammaTau);
+ double dy = 2.0 * a * a * effStressTpdt + dGammaTau *
+ (2.0 * a * alpha * dt * effStressTpdt * effStressTpdt +
+ c - 2.0 * d * d * gammaTau);
PetscLogFlops(36);
return dy;
} // _effStressDFunc
@@ -443,34 +442,35 @@
// Effective stress function that computes effective stress function
// and derivative.
void
-pylith::materials::PowerLaw3D::_effStressFuncDFunc(double x,
- void *params,
+pylith::materials::PowerLaw3D::_effStressFuncDFunc(double effStressTpdt,
+ double *params,
double *y,
double *dy)
{ // _effStressFunc
- struct effStressParams *p = (struct effStressParams *) params;
- double ae = p->ae;
- double b = p->b;
- double c = p->c;
- double d = p->d;
- double alfa = p->alfa;
- double dt = p->dt;
- double effStressT = p->effStressT;
- double powerLawExp = p->powerLawExp;
- double viscosityCoeff = p->viscosityCoeff;
- double factor1 = 1.0-alfa;
- double effStressTau = factor1 * effStressT + alfa * x;
+ double ae = params[0];
+ double b = params[1];
+ double c = params[2];
+ double d = params[3];
+ double alpha = params[4];
+ double dt = params[5];
+ double effStressT = params[6];
+ double powerLawExp = params[7];
+ double viscosityCoeff = params[8];
+ double factor1 = 1.0-alpha;
+ double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
double gammaTau = 0.5 * ((effStressTau/viscosityCoeff)^
(powerLawExp - 1.0))/viscosityCoeff;
- double dGammaTau = 0.5 * alfa * (powerLawExp - 1.0) *
+ double dGammaTau = 0.5 * alpha * (powerLawExp - 1.0) *
((effStressTau/viscosityCoeff) ^ (powerLawExp - 2.0))/
(viscosityCoeff * viscosityCoeff);
- double a = ae + alfa * dt * gammaTau;
- double *y = a * a * x * x - b + c * gammaTau - d * d * gammaTau * gammaTau;
- double *dy = 2.0 * a * a * x + dGammaTau *
- (2.0 * a * alfa * dt * x * x + c - 2.0 * d * d * gammaTau);
+ double a = ae + alpha * dt * gammaTau;
+ double *y = a * a * effStressTpdt * effStressTpdt - b +
+ c * gammaTau - d * d * gammaTau * gammaTau;
+ double *dy = 2.0 * a * a * effStressTpdt + dGammaTau *
+ (2.0 * a * alpha * dt * effStressTpdt * effStressTpdt +
+ c - 2.0 * d * d * gammaTau);
PetscLogFlops(46);
-} // _effStressFunc
+} // _effStressFuncDFunc
// ----------------------------------------------------------------------
// Compute stress tensor at location from properties as a viscoelastic
// material.
@@ -500,18 +500,100 @@
const double lambda = properties[_PowerLaw3D::pidLambda];
const double viscosityCoeff = properties[_PowerLaw3D::pidViscosityCoeff];
const double powerLawExp = properties[_PowerLaw3D::pidPowerLawExp];
+ memcpy(&_visStrainT[0], &properties[_PowerLaw3D::pidVisStrainT],
+ tensorSize * sizeof(double));
+ memcpy(&_stressT[0], &properties[_PowerLaw3D::pidStressT],
+ tensorSize * sizeof(double));
const double mu2 = 2.0 * mu;
+ const double lamPlusMu = lambda + mu;
const double bulkModulus = lambda + mu2/3.0;
+ const double youngs = mu * (3.0 * lambda + mu2)/lamPlusMu;
+ const double poissons = 0.5 * lambda/lamPlusMu;
+ const double ae = (1.0 + poissons)/youngs;
+ // Need to figure out how time integration parameter alpha is going to be
+ // specified. It should probably be specified in the problem definition and
+ // then used only by the material types that use it. For now we are setting
+ // it to 0.5, which should probably be the default value.
+ const double alpha = 0.5;
+ const double timeFac = _dt * (1.0 - alpha);
+
+ // Values for current time step
const double e11 = totalStrain[0];
const double e22 = totalStrain[1];
const double e33 = totalStrain[2];
-
const double traceStrainTpdt = e11 + e22 + e33;
const double meanStrainTpdt = traceStrainTpdt/3.0;
const double meanStressTpdt = bulkModulus * traceStrainTpdt;
+ const double devStrainTpdt[] = { _totalStrain[0] - meanStrainTpdt,
+ _totalStrain[1] - meanStrainTpdt,
+ _totalStrain[2] - meanStrainTpdt,
+ _totalStrain[3],
+ _totalStrain[4],
+ _totalStrain[5] };
+ const double strainInvar2Tpdt = 0.5 *
+ _scalarProduct(devStrainTpdt, devStrainTpdt);
+ // Values for previous time step
+ const double meanStressT = (_stressT[0] + _stressT[1] + _stressT[2])/3.0;
+ const double devStressT[] = { _stressT[0] - meanStressT,
+ _stressT[1] - meanStressT,
+ _stressT[2] - meanStressT,
+ _stressT[3],
+ _stressT[4],
+ _stressT[5] };
+ const double stressInvar2T = 0.5 * _scalarProduct(devStressT, devStressT);
+ const double effStressT = sqrt(stressInvar2T);
+
+ // Initial stress values
+ const double meanStressInitial = (_stressInitial[0] + _stressInitial[1] +
+ _stressInitial[2])/3.0;
+ const double devStressInitial[] = { _stressInitial[0] - meanStressInitial,
+ _stressInitial[1] - meanStressInitial,
+ _stressInitial[2] - meanStressInitial,
+ _stressInitial[3],
+ _stressInitial[4],
+ _stressInitial[5] };
+ const double stressInvar2Initial = 0.5 *
+ _scalarProduct(devStressInitial, devStressInitial);
+
+ // Finish defining parameters needed for root-finding algorithm.
+ const double b = strainInvar2Tpdt +
+ ae * _scalarProduct(devStrainTpdt, devStressInitial) +
+ ae * ae * stressInvar2Initial;
+ const double c = (_scalarProduct(devStrainTpdt, devStressT) +
+ ae * _scalarProduct(devStressT, devStressInitial)) * timeFac;
+ const double d = timeFac * effStressT;
+
+ // Put parameters into a vector and call root-finding algorithm.
+ // This could also be a struct.
+ const double effStressParams[] = {ae,
+ b,
+ c,
+ d,
+ alpha,
+ _dt,
+ effectiveStressT,
+ powerLawExp,
+ viscosityCoeff};
+ // I think the PETSc root-finding procedure is too involved for what we want
+ // here. I would like the function to work something like:
+ const double effStressInitialGuess = effStressT;
+ double effStressTpdt = effStressRoot(effStressInitialGuess,
+ effStressParams,
+ effStressFunc,
+ effStressFuncDFunc);
+ // I think it would be pretty easy to write a 1D root-finding algorithm that
+ // combines Newton with bisection. It would first have to bracket the root,
+ // then use Newton unless that throws the solution out of bounds or is moving
+ // too slowly.
+
+ // ********** Need to finish fixing things from here. Once I have the
+ // effective stress, I can use it to compute the current stress estimates,
+ // etc. ************
+
+
const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
// Get viscous strains
@@ -523,7 +605,7 @@
initialState,
initialStateSize);
} else {
- memcpy(&_visStrain[0], &properties[_PowerLaw3D::pidVisStrain],
+ memcpy(&_visStrainT[0], &properties[_PowerLaw3D::pidVisStrainT],
tensorSize * sizeof(double));
} // else
@@ -531,7 +613,7 @@
double devStressTpdt = 0.0;
for (int iComp=0; iComp < tensorSize; ++iComp) {
- devStressTpdt = mu2 * _visStrain[iComp];
+ devStressTpdt = mu2 * _visStrainT[iComp];
// Later I will want to put in initial stresses.
stress[iComp] = diag[iComp] * meanStressTpdt + devStressTpdt +
@@ -695,7 +777,7 @@
for (int iComp=0; iComp < _PowerLaw3D::tensorSize; ++iComp) {
properties[_PowerLaw3D::pidStrainT+iComp] = totalStrain[iComp];
- properties[_PowerLaw3D::pidVisStrain+iComp] =
+ properties[_PowerLaw3D::pidVisStrainT+iComp] =
totalStrain[iComp] - diag[iComp] * meanStrainTpdt;
} // for
PetscLogFlops(3 + 2 * _PowerLaw3D::tensorSize);
@@ -729,8 +811,8 @@
initialState,
initialStateSize);
- memcpy(&properties[_PowerLaw3D::pidVisStrain],
- &_visStrain[0],
+ memcpy(&properties[_PowerLaw3D::pidVisStrainT],
+ &_visStrainT[0],
tensorSize * sizeof(double));
memcpy(&properties[_PowerLaw3D::pidStrainT],
&totalStrain[0],
Modified: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh 2009-04-23 02:32:01 UTC (rev 14788)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh 2009-04-24 04:52:48 UTC (rev 14789)
@@ -10,44 +10,45 @@
// ----------------------------------------------------------------------
//
-/** @file libsrc/materials/MaxwellIsotropic3D.h
+/** @file libsrc/materials/PowerLaw3D.h
*
- * @brief C++ MaxwellIsotropic3D object
+ * @brief C++ PowerLaw3D object
*
- * 3-D, isotropic, linear Maxwell viscoelastic material. The
+ * 3-D, isotropic, power-law Maxwell viscoelastic material. The
* physical properties are specified using density, shear-wave speed,
- * viscosity, and compressional-wave speed. The physical properties are
- * stored internally using density, lambda, mu, which are directly
- * related to the elasticity constants used in the finite-element
- * integration. The viscosity is stored using Maxwell Time (viscosity/mu).
+ * viscosity coefficient, power-law exponent, and compressional-wave speed.
+ * The physical properties are stored internally using density, lambda, mu,
+ * which are directly related to the elasticity constants used in the
+ * finite-element integration. The viscosity information is retained as
+ * specified.
*/
-#if !defined(pylith_materials_maxwellisotropic3d_hh)
-#define pylith_materials_maxwellisotropic3d_hh
+#if !defined(pylith_materials_powerlaw3d_hh)
+#define pylith_materials_powerlaw3d_hh
#include "ElasticMaterial.hh"
/// Namespace for pylith package
namespace pylith {
namespace materials {
- class MaxwellIsotropic3D;
- class TestMaxwellIsotropic3D; // unit testing
+ class PowerLaw3D;
+ class TestPowerLaw3D; // unit testing
} // materials
} // pylith
/// 3-D, isotropic, linear Maxwell viscoelastic material.
-class pylith::materials::MaxwellIsotropic3D : public ElasticMaterial
-{ // class MaxwellIsotropic3D
- friend class TestMaxwellIsotropic3D; // unit testing
+class pylith::materials::PowerLaw3D : public ElasticMaterial
+{ // class PowerLaw3D
+ friend class TestPowerLaw3D; // unit testing
// PUBLIC METHODS /////////////////////////////////////////////////////
public :
/// Default constructor
- MaxwellIsotropic3D(void);
+ PowerLaw3D(void);
/// Destructor
- ~MaxwellIsotropic3D(void);
+ ~PowerLaw3D(void);
/** Set current time step.
*
@@ -195,7 +196,7 @@
private :
/// Member prototype for _calcStress()
- typedef void (pylith::materials::MaxwellIsotropic3D::*calcStress_fn_type)
+ typedef void (pylith::materials::PowerLaw3D::*calcStress_fn_type)
(double* const,
const int,
const double*,
@@ -207,7 +208,7 @@
const bool);
/// Member prototype for _calcElasticConsts()
- typedef void (pylith::materials::MaxwellIsotropic3D::*calcElasticConsts_fn_type)
+ typedef void (pylith::materials::PowerLaw3D::*calcElasticConsts_fn_type)
(double* const,
const int,
const double*,
@@ -218,7 +219,7 @@
const int);
/// Member prototype for _updateProperties()
- typedef void (pylith::materials::MaxwellIsotropic3D::*updateProperties_fn_type)
+ typedef void (pylith::materials::PowerLaw3D::*updateProperties_fn_type)
(double* const,
const int,
const double*,
@@ -367,10 +368,10 @@
private :
/// Not implemented
- MaxwellIsotropic3D(const MaxwellIsotropic3D& m);
+ PowerLaw3D(const PowerLaw3D& m);
/// Not implemented
- const MaxwellIsotropic3D& operator=(const MaxwellIsotropic3D& m);
+ const PowerLaw3D& operator=(const PowerLaw3D& m);
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
@@ -387,11 +388,11 @@
/// Method to use for _updateProperties().
updateProperties_fn_type _updatePropertiesFn;
-}; // class MaxwellIsotropic3D
+}; // class PowerLaw3D
-#include "MaxwellIsotropic3D.icc" // inline methods
+#include "PowerLaw3D.icc" // inline methods
-#endif // pylith_materials_maxwellisotropic3d_hh
+#endif // pylith_materials_powerlaw3d_hh
// End of file
Modified: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.icc 2009-04-23 02:32:01 UTC (rev 14788)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.icc 2009-04-24 04:52:48 UTC (rev 14789)
@@ -10,8 +10,8 @@
// ----------------------------------------------------------------------
//
-#if !defined(pylith_materials_maxwellisotropic3d_hh)
-#error "MaxwellIsotropic3D.icc can only be included from MaxwellIsotropic3D.hh"
+#if !defined(pylith_materials_powerlaw3d_hh)
+#error "PowerLaw3D.icc can only be included from PowerLaw3D.hh"
#endif
#include <assert.h> // USES assert()
@@ -20,45 +20,51 @@
// Set current time step.
inline
void
-pylith::materials::MaxwellIsotropic3D::timeStep(const double dt) {
- // Jacobian needs to be reformed if the time step size changes.
- if (_dt > 0.0 && dt != _dt)
- _needNewJacobian = true;
+pylith::materials::PowerLaw3D::timeStep(const double dt) {
+ // Not sure what to do here. If we are using full Newton the Jacobian will
+ // always need reforming, but SNES may opt not to reform it sometimes.
+ _needNewJacobian = true;
_dt = dt;
} // timeStep
// Set whether elastic or inelastic constitutive relations are used.
inline
void
-pylith::materials::MaxwellIsotropic3D::useElasticBehavior(const bool flag) {
+pylith::materials::PowerLaw3D::useElasticBehavior(const bool flag,
+ const int iteration) {
if (flag) {
_calcStressFn =
- &pylith::materials::MaxwellIsotropic3D::_calcStressElastic;
+ &pylith::materials::PowerLaw3D::_calcStressElastic;
_calcElasticConstsFn =
- &pylith::materials::MaxwellIsotropic3D::_calcElasticConstsElastic;
+ &pylith::materials::PowerLaw3D::_calcElasticConstsElastic;
_updatePropertiesFn =
- &pylith::materials::MaxwellIsotropic3D::_updatePropertiesElastic;
+ &pylith::materials::PowerLaw3D::_updatePropertiesElastic;
} else {
_calcStressFn =
- &pylith::materials::MaxwellIsotropic3D::_calcStressViscoelastic;
- _calcElasticConstsFn =
- &pylith::materials::MaxwellIsotropic3D::_calcElasticConstsViscoelastic;
+ &pylith::materials::PowerLaw3D::_calcStressViscoelastic;
+ if (iteration == 0) {
+ _calcElasticConstsFn =
+ &pylith::materials::PowerLaw3D::_calcElasticConstsViscoelasticInitial;
+ } else {
+ _calcElasticConstsFn =
+ &pylith::materials::PowerLaw3D::_calcElasticConstsViscoelastic;
+ }
_updatePropertiesFn =
- &pylith::materials::MaxwellIsotropic3D::_updatePropertiesViscoelastic;
+ &pylith::materials::PowerLaw3D::_updatePropertiesViscoelastic;
} // if/else
} // useElasticBehavior
// Get flag indicating whether material implements an empty
inline
bool
-pylith::materials::MaxwellIsotropic3D::usesUpdateProperties(void) const {
+pylith::materials::PowerLaw3D::usesUpdateProperties(void) const {
return true;
} // usesUpdateProperties
// Compute stress tensor from parameters.
inline
void
-pylith::materials::MaxwellIsotropic3D::_calcStress(double* const stress,
+pylith::materials::PowerLaw3D::_calcStress(double* const stress,
const int stressSize,
const double* parameters,
const int numParams,
@@ -78,7 +84,7 @@
// Compute derivatives of elasticity matrix from parameters.
inline
void
-pylith::materials::MaxwellIsotropic3D::_calcElasticConsts(
+pylith::materials::PowerLaw3D::_calcElasticConsts(
double* const elasticConsts,
const int numElasticConsts,
const double* parameters,
@@ -97,7 +103,7 @@
// Update state variables after solve.
inline
void
-pylith::materials::MaxwellIsotropic3D::_updateProperties(double* const parameters,
+pylith::materials::PowerLaw3D::_updateProperties(double* const parameters,
const int numParams,
const double* totalStrain,
const int strainSize,
@@ -109,4 +115,19 @@
initialState, initialStateSize);
} // _updateProperties
+// Compute scalar product, assuming vector form of a tensor.
+inline
+double
+pylith::materials::PowerLaw3D::_scalarProduct(const double* tensor1,
+ const double* tensor2)
+{
+ const double scalarProduct = tensor1[0] * tensor2[0] +
+ tensor1[1] * tensor2[1] +
+ tensor1[2] * tensor2[2] +
+ 2.0 * (tensor1[3] * tensor2[3] +
+ tensor1[4] * tensor2[4] +
+ tensor1[5] * tensor2[5]);
+ return scalarProduct;
+} // _scalarProduct
+
// End of file
More information about the CIG-COMMITS
mailing list