C++ MDArray::rank方法代码示例

void IntrepidSideCell<MDArray>::mapToCellPhysicalFrame( 
    const MDArray& parametric_coords, MDArray& physical_coords )
{
    DTK_REQUIRE( 2 == parametric_coords.rank() );
    DTK_REQUIRE( 3 == physical_coords.rank() );
    DTK_REQUIRE( parametric_coords.dimension(1) ==
		   Teuchos::as<int>(this->d_topology.getDimension()) );
    DTK_REQUIRE( physical_coords.dimension(0) ==
		   this->d_cell_node_coords.dimension(0) );
    DTK_REQUIRE( physical_coords.dimension(1) ==
		   parametric_coords.dimension(0) );
    DTK_REQUIRE( physical_coords.dimension(2) ==
		   Teuchos::as<int>(d_parent_topology.getDimension()) );

    MDArray mapped_coords( parametric_coords.dimension(0),
			   d_parent_topology.getDimension() );

    Intrepid::CellTools<Scalar>::mapToReferenceSubcell(
	mapped_coords, parametric_coords, this->d_topology.getDimension(),
	d_side_id, d_parent_topology );

    Intrepid::CellTools<Scalar>::mapToPhysicalFrame( 
	physical_coords, mapped_coords, 
	this->d_cell_node_coords, d_parent_topology );
}

    void StringFunction::operator()(MDArray& input_phy_points, MDArray& output_values,
                                    const stk::mesh::Bucket& bucket, const MDArray& parametric_coordinates, double time_value_optional) 
    {
      PRINT("tmp srk StringFunction::operator(bucket) getName()= " << getName() << " input_phy_points= " << input_phy_points << " output_values= " << output_values);

      VERIFY_OP(input_phy_points.rank(), ==, 3, "StringFunction::operator() must pass in input_phy_points(numCells, numPointsPerCell, spaceDim)");
      int nPoints = input_phy_points.dimension(1);
      int nSpaceDim = input_phy_points.dimension(2);
      int nOutDim = output_values.dimension(2);
      MDArray input_phy_points_one(1, nPoints, nSpaceDim);
      MDArray output_values_one   (1, nPoints, nOutDim);
      const unsigned num_elements_in_bucket = bucket.size();
      for (unsigned iElement = 0; iElement < num_elements_in_bucket; iElement++)
        {
          stk::mesh::Entity& element = bucket[iElement];
          for (int iPoint = 0; iPoint<nPoints; iPoint++)
            {
              for (int iSpaceDim=0; iSpaceDim < nSpaceDim; iSpaceDim++)
                {
                  input_phy_points_one(0, iPoint, iSpaceDim) = input_phy_points(iElement, iPoint, iSpaceDim);
                }
            }
          (*this)(input_phy_points_one, output_values_one, element, parametric_coordinates, time_value_optional);
          for (int iPoint = 0; iPoint<nPoints; iPoint++)
            {
              for (int iDOF = 0; iDOF < nOutDim; iDOF++)
                {
                  output_values(iElement, iPoint, iDOF) = output_values_one(0, iPoint, iDOF);
                }
            }
        }
    }

void NewtonSolver<NonlinearProblem>::solve( MDArray& u, 
					    NonlinearProblem& problem,
					    const double tolerance,
					    const int max_iters )
{
    DTK_REQUIRE( 2 == u.rank() );

    // Allocate nonlinear residual, Jacobian, Newton update, and work arrays.
    int d0 = u.dimension(0);
    int d1 = u.dimension(1);
    MDArray F( d0, d1 );
    MDArray J( d0, d1, d1 );
    MDArray J_inv( d0, d1, d1 );
    MDArray update( d0, d1 );
    MDArray u_old = u;
    MDArray conv_check( d0 );

    // Compute the initial state.
    NPT::updateState( problem, u );

    // Computen the initial nonlinear residual and scale by -1 to get -F(u).
    NPT::evaluateResidual( problem, u, F );
    Intrepid::RealSpaceTools<Scalar>::scale( 
	F, -Teuchos::ScalarTraits<Scalar>::one() );

    // Compute the initial Jacobian.
    NPT::evaluateJacobian( problem, u, J );

    // Check for degeneracy of the Jacobian. If it is degenerate then the
    // problem is ill conditioned and return very large numbers in the state
    // vector that correspond to no solution.
    MDArray det( 1 );
    Intrepid::RealSpaceTools<Scalar>::det( det, J );
    if ( std::abs(det(0)) < tolerance )
    {
	for ( int m = 0; m < d0; ++m )
	{
	    for ( int n = 0; n < d1; ++n )
	    {
		u(m,n) = std::numeric_limits<Scalar>::max();
	    }
	}
	return;
    }

    // Nonlinear solve.
    for ( int k = 0; k < max_iters; ++k )
    {
	// Solve the linear model, delta_u = J^-1 * -F(u).
	Intrepid::RealSpaceTools<Scalar>::inverse( J_inv, J );
	Intrepid::RealSpaceTools<Scalar>::matvec( update, J_inv, F );

	// Update the solution, u += delta_u.
	Intrepid::RealSpaceTools<Scalar>::add( u, update );

	// Check for convergence.
	Intrepid::RealSpaceTools<Scalar>::subtract( u_old, u );
	Intrepid::RealSpaceTools<Scalar>::vectorNorm( 
	    conv_check, u_old, Intrepid::NORM_TWO );
	if ( tolerance > conv_check(0) )
	{
	    break;
	}

	// Reset for the next iteration.
	u_old = u;

	// Update any state-dependent data from the last iteration using the
	// new solution vector.
	NPT::updateState( problem, u );

	// Compute the nonlinear residual and scale by -1 to get -F(u).
	NPT::evaluateResidual( problem, u, F );
	Intrepid::RealSpaceTools<Scalar>::scale( 
	    F, -Teuchos::ScalarTraits<Scalar>::one() );

	// Compute the Jacobian.
	NPT::evaluateJacobian( problem, u, J );
    }

    // Check for convergence.
    DTK_ENSURE( tolerance > conv_check(0) );
}

    /** Evaluate the function at this input point (or points) returning value(s) in output_field_values
     *
     *   In the following, the arrays are dimensioned using the notation (from Intrepid's doc):
     *
     *   [C]         - num. integration domains (cells/elements)
     *   [F]         - num. Intrepid "fields" (number of bases within an element == num. nodes typically)
     *   [P]         - num. integration (or interpolation) points within the element
     *   [D]         - spatial dimension
     *   [D1], [D2]  - spatial dimension
     *
     *   Locally, we introduce this notation:
     *
     *   [DOF]       - number of degrees-of-freedom per node of the interpolated stk Field.  For example, a vector field in 3D has [DOF] = 3
     *
     *  Dimensions of input_phy_points are required to be either ([D]) or ([P],[D])
     *  Dimensions of output_field_values are required to be ([DOF]) or ([P],[DOF]) respectively
     *
     *  [R] is used for the rank of MDArray's
     */
    void FieldFunction::operator()(MDArray& input_phy_points, MDArray& output_field_values, double time)
    {
      EXCEPTWATCH;
      argsAreValid(input_phy_points, output_field_values);

      m_found_on_local_owned_part = false;

      //// single point only (for now)
      unsigned found_it = 0;

      int D_ = last_dimension(input_phy_points);
      MDArray found_parametric_coordinates_one(1, D_);
      setup_searcher(D_);

      MDArray output_field_values_local = output_field_values;
      int R_output = output_field_values.rank();

      int R_input = input_phy_points.rank();
      int P_ = (R_input == 1 ? 1 : input_phy_points.dimension(R_input-2));

      // FIXME for tensor valued fields
      int DOF_ = last_dimension(output_field_values_local);

      MDArray input_phy_points_one(1,D_);
      MDArray output_field_values_one(1,DOF_);

      int C_ = 1;
      if (R_input == 3)
        {
          C_ = input_phy_points.dimension(0);
        }
      for (int iC = 0; iC < C_; iC++)
        {
          for (int iP = 0; iP < P_; iP++)
            {
              for (int iD = 0; iD < D_; iD++)
                {
                  switch(R_input)
                    {
                    case 1: input_phy_points_one(0, iD) = input_phy_points(iD); break;
                    case 2: input_phy_points_one(0, iD) = input_phy_points(iP, iD); break;
                    case 3: input_phy_points_one(0, iD) = input_phy_points(iC, iP, iD); break;
                    default: VERIFY_1("bad rank");
                    }
                }

              const stk_classic::mesh::Entity *found_element = 0;
              {
                EXCEPTWATCH;
                //if (m_searchType==STK_SEARCH) std::cout << "find" << std::endl;
                found_element = m_searcher->findElement(input_phy_points_one, found_parametric_coordinates_one, found_it, m_cachedElement);
                //if (m_searchType==STK_SEARCH)                std::cout << "find..done found_it=" << found_it << std::endl;
              }

              // if found element on the local owned part, evaluate
              if (found_it)
                {
                  m_found_on_local_owned_part = true;
                  if (( EXTRA_PRINT) && m_searchType==STK_SEARCH)
                    std::cout << "FieldFunction::operator() found element # = " << found_element->identifier() << std::endl;

                  (*this)(input_phy_points_one, output_field_values_one, *found_element, found_parametric_coordinates_one);

                  for (int iDOF = 0; iDOF < DOF_; iDOF++)
                    {
                      switch (R_output)
                        {
                        case 1: output_field_values_local( iDOF)        = output_field_values_one(0, iDOF); break;
                        case 2: output_field_values_local(iP, iDOF)     = output_field_values_one(0, iDOF); break;
                        case 3: output_field_values_local(iC, iP, iDOF) = output_field_values_one(0, iDOF); break;
                        default: VERIFY_1("bad rank");
                        }
                    }
                }
              else
                {
                  if (!m_parallelEval)
                    {
                      std::cout << "P[" << Util::get_rank() << "] FieldFunction::operator() found_it = " << found_it << " points= "
                                << input_phy_points_one
                                << std::endl;
//.........这里部分代码省略.........