自学教程：C++ tensor函数代码示例

void InertiaSensorFilter::covOfSigmaPoints(){  cov = tensor(sigmaPoints[0] - x) +        (tensor(sigmaPoints[1] - x) + tensor(sigmaPoints[2] - x)) +        (tensor(sigmaPoints[3] - x) + tensor(sigmaPoints[4] - x));  cov *= 0.5f;}

static void test_static_index_list(){    Tensor<float, 4> tensor(2,3,5,7);    tensor.setRandom();    constexpr auto reduction_axis = make_index_list(0, 1, 2);    VERIFY_IS_EQUAL(internal::array_get<0>(reduction_axis), 0);    VERIFY_IS_EQUAL(internal::array_get<1>(reduction_axis), 1);    VERIFY_IS_EQUAL(internal::array_get<2>(reduction_axis), 2);    VERIFY_IS_EQUAL(static_cast<DenseIndex>(reduction_axis[0]), 0);    VERIFY_IS_EQUAL(static_cast<DenseIndex>(reduction_axis[1]), 1);    VERIFY_IS_EQUAL(static_cast<DenseIndex>(reduction_axis[2]), 2);    EIGEN_STATIC_ASSERT((internal::array_get<0>(reduction_axis) == 0), YOU_MADE_A_PROGRAMMING_MISTAKE);    EIGEN_STATIC_ASSERT((internal::array_get<1>(reduction_axis) == 1), YOU_MADE_A_PROGRAMMING_MISTAKE);    EIGEN_STATIC_ASSERT((internal::array_get<2>(reduction_axis) == 2), YOU_MADE_A_PROGRAMMING_MISTAKE);    Tensor<float, 1> result = tensor.sum(reduction_axis);    for (int i = 0; i < result.size(); ++i) {        float expected = 0.0f;        for (int j = 0; j < 2; ++j) {            for (int k = 0; k < 3; ++k) {                for (int l = 0; l < 5; ++l) {                    expected += tensor(j,k,l,i);                }            }        }        VERIFY_IS_APPROX(result(i), expected);    }}

		/**		Returns nullptr if settings incorrect.		*/		std::unique_ptr<Configurator> FractureConfiguratorFactory::configurator(const Settings& s) const		{			std::unique_ptr<Configurator> pConf(nullptr);			const string c = s.json["configuration"].get<string>();			if (c == string("uniform")) {				const double am = s.json["mechanical aperture"].get<double>();				if (s.json["hydraulic aperture"].size() == 9) // tensor					pConf.reset(new UniformFractureConfigurator(tensor("hydraulic aperture", s),																tensor("permeability", s), 																tensor("conductivity", s),																am));				else // scalar					pConf.reset(new UniformFractureConfigurator(s.json["hydraulic aperture"].get<double>(), am));			}			else if (c == string("regional uniform")) {				const auto ams = s.json["mechanical aperture"].get<vector<double>>();				const auto frnames = s.json["fracture regions"].get<vector<string>>();				if (s.json.count("permeability")) { // tensor					vector<string> props = { "hydraulic aperture" , "permeability", "conductivity" };					vector<vector<TensorVariable<3>>> vals(props.size(), vector<TensorVariable<3>>(ams.size(), TensorVariable<3>()));					for (size_t i(0); i < props.size(); ++i) {						auto jvals = s.json[props[i].c_str()].get<vector<vector<double>>>();						for (size_t j(0); j < ams.size(); ++j)							vals[i][j] = tensor(jvals.at(j));					}					pConf.reset(new RegionalUniformFractureConfigurator(vals.at(0), vals.at(1), vals.at(2), ams, frnames));				}				else { // scalar					const auto ahs = s.json["hydraulic aperture"].get<vector<double>>();					pConf.reset(new RegionalUniformFractureConfigurator(ahs, ams, frnames));				}			}			return pConf;		}

static void test_dynamic_index_list(){    Tensor<float, 4> tensor(2,3,5,7);    tensor.setRandom();    int dim1 = 2;    int dim2 = 1;    int dim3 = 0;    auto reduction_axis = make_index_list(dim1, dim2, dim3);    VERIFY_IS_EQUAL(internal::array_get<0>(reduction_axis), 2);    VERIFY_IS_EQUAL(internal::array_get<1>(reduction_axis), 1);    VERIFY_IS_EQUAL(internal::array_get<2>(reduction_axis), 0);    VERIFY_IS_EQUAL(static_cast<DenseIndex>(reduction_axis[0]), 2);    VERIFY_IS_EQUAL(static_cast<DenseIndex>(reduction_axis[1]), 1);    VERIFY_IS_EQUAL(static_cast<DenseIndex>(reduction_axis[2]), 0);    Tensor<float, 1> result = tensor.sum(reduction_axis);    for (int i = 0; i < result.size(); ++i) {        float expected = 0.0f;        for (int j = 0; j < 2; ++j) {            for (int k = 0; k < 3; ++k) {                for (int l = 0; l < 5; ++l) {                    expected += tensor(j,k,l,i);                }            }        }        VERIFY_IS_APPROX(result(i), expected);    }}

void InertialDataFilter::covOfSigmaPoints(){  cov = tensor(sigmaPoints[0] - mean) +    tensor(sigmaPoints[1] - mean) +    tensor(sigmaPoints[2] - mean) +    tensor(sigmaPoints[3] - mean) +    tensor(sigmaPoints[4] - mean);  cov *= 0.5f;}

Matrix2x2<> AngleEstimator::covOfSigmaPZgyroAndSigmaPX(const Vector2<>& mean) const{  Matrix2x2<> result(tensor(sigmaPZgyro[1] - mean, l.c[0]));  result += tensor(sigmaPZgyro[2] - mean, l.c[1]);  result += tensor(sigmaPZgyro[3] - mean, -l.c[0]);  result += tensor(sigmaPZgyro[4] - mean, -l.c[1]);  result /= 2.;  return result;}

void Foam::coordinateRotationOFext::calcTransform(    const vector& axis1,    const vector& axis2,    const axisOrder& order){    vector a = axis1 / mag(axis1);    vector b = axis2;    // Absorb minor non-orthogonality into axis2    b = b - (b & a)*a;    if (mag(b) < SMALL)    {        FatalErrorIn("coordinateRotation::calcTransform()")            << "axis1, axis2 appear co-linear: "            << axis1 << ", " << axis2 << endl            << abort(FatalError);    }    b = b / mag(b);    vector c = a ^ b;    // the global -> local transformation    tensor Rtr;    switch (order)    {        case e1e2:            Rtr = tensor(a, b, c);            break;        case e2e3:            Rtr = tensor(c, a, b);            break;        case e3e1:            Rtr = tensor(b, c, a);            break;        default:            FatalErrorIn("coordinateRotation::calcTransform()")                << "programmer error" << endl                << abort(FatalError);            // To satisfy compiler warnings            Rtr = tensor::zero;            break;    }    // the local -> global transformation    tensor::operator=( Rtr.T() );}

bool MxQuadric3::optimize(Vec3& v, const Vec3& v1,			  const Vec3& v2, const Vec3& v3) const{    Vec3 d13 = v1 - v3;    Vec3 d23 = v2 - v3;    Mat3 A = tensor();    Vec3 B = vector();    Vec3 Ad13 = A*d13;    Vec3 Ad23 = A*d23;    Vec3 Av3  = A*v3;    double d13_d23 = (d13*Ad23) + (d23*Ad13);    double v3_d13  = (d13*Av3) + (v3*Ad13);    double v3_d23  = (d23*Av3) + (v3*Ad23);    double d23Ad23 = d23*Ad23;    double d13Ad13 = d13*Ad13;    double denom = d13Ad13*d23Ad23 - 2*d13_d23;    if( FEQ(denom, 0.0, 1e-12) )	return false;    double a = ( d23Ad23*(2*(B*d13) + v3_d13) -		  d13_d23*(2*(B*d23) + v3_d23) ) / -denom;    double b = ( d13Ad13*(2*(B*d23) + v3_d23) -		  d13_d23*(2*(B*d13) + v3_d13) ) / -denom;    if( a<0.0 ) a=0.0; else if( a>1.0 ) a=1.0;    if( b<0.0 ) b=0.0; else if( b>1.0 ) b=1.0;    v = a*d13 + b*d23 + v3;    return true;}

EPState::EPState( )//ZC05/2004: Eo(30000.0), E_Young(30000.0), nu_Poisson(0.3), rho_mass_density(0.0),//ZC05/2004  Converged(false),//ZC05/2004  Elasticflag(0),Ev(0.0),nuhv(0.0),Ghv(0.0),//ZC05/2004  eo(0.85), ec(0.80), Lambda(0.025), po(100.0), e(0.85), psi(0.05), a(0.5): Converged(false), e(0.85), psi(0.05) //ZC05/2004{      //Eo               = 30000.0;      //E_Young          = 30000.0;      //nu_Poisson       = 0.3;      //rho_mass_density = 0.0;      Eep = tensor( 4, def_dim_4, 0.0 );      integratorFlag = 0;//ForwardEuler assumed      Delta_lambda = 0.0;      NScalarVar = MaxNScalarVar;      for (int i =0; i < NScalarVar; i++) {      	 ScalarVar_commit[i] = 0.0;         ScalarVar[i] = 0.0;	 ScalarVar_init[i] = 0.0;      }            NTensorVar = MaxNTensorVar;      for (int j =0; j < NTensorVar; j++) {         TensorVar[j] = stresstensor(0.0);         TensorVar_commit[j] = stresstensor(0.0);         TensorVar_init[j] = stresstensor(0.0);      }      //Converged = false;}

Tensor variable_data_factory(const Type& type, PyObject* args, PyObject* kwargs) {  static PythonArgParser parser({    "new(Tensor other, *, bool requires_grad=False)",    "new(PyObject* data, *, int64_t device=-1, bool requires_grad=False)",  });  PyObject* parsed_args[3];  auto r = parser.parse(args, kwargs, parsed_args);  if (r.idx == 0) {    return set_requires_grad(new_with_tensor_copy(type, r.tensor(0)), r.toBool(1));    return set_requires_grad(new_with_tensor(type, r.tensor(0)), r.toBool(1));  } else if (r.idx == 1) {    return set_requires_grad(new_from_data(type, r.toInt64(1), r.pyobject(0)), r.toBool(2));  }  throw std::runtime_error("variable(): invalid arguments");}

Tensor tensor_new(const Type& type, PyObject* args, PyObject* kwargs) {  static PythonArgParser parser({    "new(*, int64_t device=-1)",    "new(IntList size, *, int64_t device=-1)",    "new(Storage storage)",    "new(*, int64_t cdata)|hidden",    "new(Tensor other)",    "new(PyObject* data, *, int64_t device=-1)",  });  PyObject* parsed_args[2];  auto r = parser.parse(args, kwargs, parsed_args);  if (r.idx == 0) {    AutoGPU auto_gpu(r.toInt64(0));    return type.tensor();  } else if (r.idx == 1) {    PyObject* arg = parsed_args[0];    if (!THPSize_Check(arg) && PyTuple_GET_SIZE(args) >= 1 && arg == PyTuple_GET_ITEM(args, 0)) {      // new(sequence) binds to this signature but should be treated differently      // unless the sequences is a torch.Size      return new_from_sequence(type, r.toInt64(1), r.pyobject(0));    }    return new_with_sizes(type, r.toInt64(1), r.intlist(0));  } else if (r.idx == 2) {    return new_with_storage(type, *r.storage(0));  } else if (r.idx == 3) {    auto cdata = reinterpret_cast<void*>(r.toInt64(0));    return type.unsafeTensorFromTH(cdata, true);  } else if (r.idx == 4) {    return new_with_tensor(type, r.tensor(0));  } else if (r.idx == 5) {    return new_from_sequence(type, r.toInt64(1), r.pyobject(0));  }  throw std::runtime_error("new(): invalid arguments");}

void tree1(int k, int i, int matrix[7][4], int E[7][4][4], int pivots[7][4],              int rge[7][4], int u[7], int s[7], int zero[3]){    //lookup table for indices  static int first[3][2] = {{0,1}, {0,2}, {2,2}};  static int index[17][2] =  {    {0,1}, {1,0}, {0,0}, {1,1},    {2,0}, {3,0}, {2,1}, {3,1},    {0,2}, {1,2}, {0,3}, {1,3},    {2,2}, {3,2}, {2,3}, {3,3},    {4,4}  };    //set u,s indices  if (k == 0) {    u[k] = first[i][0];    s[k] = first[i][1];  }  else {    u[k] = index[i][0];    s[k] = index[i][1];  }    //handle zero case  zero[k] = (i == 16);  tensor(matrix[k], u[k], s[k]);  //fill raw matrix  RGE_Solver(k, matrix[k], E, rge[k], pivots);  //Running Gaussian Elimination}

void btCompoundShape::calculatePrincipalAxisTransform(btScalar* masses, btTransform& principal, btVector3& inertia) const{	int n = m_children.size();	btScalar totalMass = 0;	btVector3 center(0, 0, 0);	int k;	for (k = 0; k < n; k++)	{		btAssert(masses[k]>0);		center += m_children[k].m_transform.getOrigin() * masses[k];		totalMass += masses[k];	}	btAssert(totalMass>0);	center /= totalMass;	principal.setOrigin(center);	btMatrix3x3 tensor(0, 0, 0, 0, 0, 0, 0, 0, 0);	for ( k = 0; k < n; k++)	{		btVector3 i;		m_children[k].m_childShape->calculateLocalInertia(masses[k], i);		const btTransform& t = m_children[k].m_transform;		btVector3 o = t.getOrigin() - center;		//compute inertia tensor in coordinate system of compound shape		btMatrix3x3 j = t.getBasis().transpose();		j[0] *= i[0];		j[1] *= i[1];		j[2] *= i[2];		j = t.getBasis() * j;		//add inertia tensor		tensor[0] += j[0];		tensor[1] += j[1];		tensor[2] += j[2];		//compute inertia tensor of pointmass at o		btScalar o2 = o.length2();		j[0].setValue(o2, 0, 0);		j[1].setValue(0, o2, 0);		j[2].setValue(0, 0, o2);		j[0] += o * -o.x(); 		j[1] += o * -o.y(); 		j[2] += o * -o.z();		//add inertia tensor of pointmass		tensor[0] += masses[k] * j[0];		tensor[1] += masses[k] * j[1];		tensor[2] += masses[k] * j[2];	}	tensor.diagonalize(principal.getBasis(), btScalar(0.00001), 20);	inertia.setValue(tensor[0][0], tensor[1][1], tensor[2][2]);}

voidColumnMajorMatrixTest::fillMatrix(){  TensorValue<Real> tensor( 0, 0, 0, 0, 0, 0, 0, 0, 0 );  ColumnMajorMatrix & a_ref = *a;  a_ref.fill( tensor );  CPPUNIT_ASSERT( tensor(0, 0) == a_ref(0, 0) );  CPPUNIT_ASSERT( tensor(1, 0) == a_ref(1, 0) );  CPPUNIT_ASSERT( tensor(2, 0) == a_ref(2, 0) );  CPPUNIT_ASSERT( tensor(0, 1) == a_ref(0, 1) );  CPPUNIT_ASSERT( tensor(1, 1) == a_ref(1, 1) );  CPPUNIT_ASSERT( tensor(2, 1) == a_ref(2, 1) );  CPPUNIT_ASSERT( tensor(0, 2) == a_ref(0, 2) );  CPPUNIT_ASSERT( tensor(1, 2) == a_ref(1, 2) );  CPPUNIT_ASSERT( tensor(2, 2) == a_ref(2, 2) );}

Foam::energyDissipationRate::energyDissipationRate(    const word& name,    const objectRegistry& obr,    const dictionary& dict,    const bool loadFromFiles):    name_(name),    fieldName_(dict.lookup("fieldName")),    obr_(obr),    active_(true){    // Check if the available mesh is an fvMesh, otherwise deactivate    if (!isA<fvMesh>(obr_))    {        active_ = false;        WarningIn        (            "energyDissipationRate::energyDissipationRate"            "("                "const word&, "                "const objectRegistry&, "                "const dictionary&, "                "const bool"            ")"        )   << "No fvMesh available, deactivating." << nl            << endl;    }    read(dict);    if (active_)    {        const fvMesh& mesh = refCast<const fvMesh>(obr_);                volTensorField* energyDissipationRateTensorPtr        (            new volTensorField            (                IOobject                (                    fieldName_,                    mesh.time().timeName(),                    mesh,                    IOobject::NO_READ,                    IOobject::NO_WRITE                ),                mesh,                dimensionedTensor("0", sqr(dimLength)/pow(dimTime,3), tensor(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0))            )        );        mesh.objectRegistry::store(energyDissipationRateTensorPtr);    }}

Foam::tensor Foam::molecule::rotationTensorZ(scalar phi) const{    return tensor    (        Foam::cos(phi), -Foam::sin(phi), 0,        Foam::sin(phi), Foam::cos(phi), 0,        0, 0, 1    );}

bool MxQuadric3::optimize(Vec3& v) const{    Mat3 Ainv;    double det = invert(Ainv, tensor());    if( FEQ(det, 0.0, 1e-12) )	return false;    v = -(Ainv*vector());    return true;}

voidColumnMajorMatrixTest::tensorAssignOperator(){  ColumnMajorMatrix test(3, 3);  TensorValue<Real> tensor( 1, 4, 7,                            2, 5, 8,                            3, 6, 9 );  test = tensor;  CPPUNIT_ASSERT( test == *a );  CPPUNIT_ASSERT( test.numEntries() == 9 );}

void Foam::RBD::joints::Ryxz::jcalc(    joint::XSvc& J,    const scalarField& q,    const scalarField& w,    const scalarField& qDot) const{    vector qj(q.block<vector>(qIndex_));    scalar s0 = sin(qj.x());    scalar c0 = cos(qj.x());    scalar s1 = sin(qj.y());    scalar c1 = cos(qj.y());    scalar s2 = sin(qj.z());    scalar c2 = cos(qj.z());    J.X.E() = tensor    (        c2*c0 + s2*s1*s0, s2*c1,           -c2*s0 + s2*s1*c0,       -s2*c0 + c2*s1*s0, c2*c1,            s2*s0 + c2*s1*c0,        c1*s0,           -s1,               c1*c0    );    J.X.r() = Zero;    J.S = Zero;    J.S.xx() = s2*c1;    J.S.xy() = c2;    J.S.yx() = c2*c1;    J.S.yy() = -s2;    J.S.zx() = -s1;    J.S.zz() = 1;    vector qDotj(qDot.block<vector>(qIndex_));    J.v = J.S & qDotj;    J.c = spatialVector    (        c2*c1*qDotj.z()*qDotj.x()      - s2*s1*qDotj.y()*qDotj.x()      - s2*qDotj.z()*qDotj.y(),       -s2*c1*qDotj.z()*qDotj.x()      - c2*s1*qDotj.y()*qDotj.x()      - c2*qDotj.z()*qDotj.y(),       -c1*qDotj.y()*qDotj.x(),        0,        0,        0    );}

void tree3(int k, int i, int VT[7][4], int UT[7][4], int rge[7][4],                int E[7][4][4], int pivots[7][4], int inZero[3], int inU[7],                int inS[7], int inV[7], int outU[7], int outV[7], int outS[7],                int outT[7]){  static int alt[96][3] =  { //u:0, v:1, t:2, s = 4    {0,0,0}, {0,0,1}, {0,0,2}, {0,0,3}, {0,1,0}, {0,1,1}, {0,1,2}, {0,1,3},    {0,2,0}, {0,2,1}, {0,2,2}, {0,2,3}, {0,3,0}, {0,3,1}, {0,3,2}, {0,3,3},    {0,4,0}, {0,4,1}, {0,4,2}, {0,4,3}, {1,0,0}, {1,0,1}, {1,0,2}, {1,0,3},    {1,1,0}, {1,1,1}, {1,1,2}, {1,1,3}, {1,2,0}, {1,2,1}, {1,2,2}, {1,2,3},    {1,3,0}, {1,3,1}, {1,3,2}, {1,3,3}, {1,4,0}, {1,4,1}, {1,4,2}, {1,4,3},    {2,0,0}, {2,0,1}, {2,0,2}, {2,0,3}, {2,1,0}, {2,1,1}, {2,1,2}, {2,1,3},    {2,2,0}, {2,2,1}, {2,2,2}, {2,2,3}, {2,3,0}, {2,3,1}, {2,3,2}, {2,3,3},    {2,4,0}, {2,4,1}, {2,4,2}, {2,4,3}, {3,0,0}, {3,0,1}, {3,0,2}, {3,0,3},    {3,1,0}, {3,1,1}, {3,1,2}, {3,1,3}, {3,2,0}, {3,2,1}, {3,2,2}, {3,2,3},    {3,3,0}, {3,3,1}, {3,3,2}, {3,3,3}, {3,4,0}, {3,4,1}, {3,4,2}, {3,4,3},    {4,0,0}, {4,0,1}, {4,0,2}, {4,0,3}, {4,1,0}, {4,1,1}, {4,1,2}, {4,1,3},    {4,2,0}, {4,2,1}, {4,2,2}, {4,2,3}, {4,3,0}, {4,3,1}, {4,3,2}, {4,3,3}  };    if (inZero[k] ==  1) {    outU[k] = alt[i][0];    outV[k] = alt[i][1];    outS[k] = 4;    outT[k] = alt[i][2];  }  else {    outU[k] = inU[k];    outV[k] = inV[k];    outS[k] = inS[k];    outT[k] = i;  }  tensor(VT[k], outV[k], outT[k]);  tensor(UT[k], outU[k], outT[k]);  RGE_Solver(k, VT[k], E, rge[k], pivots);  }

void Foam::barycentricParticle::tetTransform(    vector& centre,    tensor& A) const{    vector base, vertex1, vertex2;    tetGeometry(centre, base, vertex1, vertex2);    A = tensor    (        base - centre,        vertex1 - centre,        vertex2 - centre    ).T();}

voidOutputTestMaterial::computeQpProperties(){  Point p = _q_point[_qp];  Real v = std::floor(_variable[_qp] + 0.5);  Real x = std::ceil(p(0));  Real y = std::ceil(p(1));  _real_property[_qp] = v + y + x + _factor;  RealVectorValue vec(v + x, v + y);  _vector_property[_qp] = vec;  RealTensorValue tensor(v, x * y, 0, -x * y, y * y);  _tensor_property[_qp] = tensor;}