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

    ///  Dimensions of parametric_coordinates and input_phy_points are required to be ([P],[D])
    /// output_values: ([P], [DOF])
    void StringFunction::operator()(MDArray& input_phy_points, MDArray& output_values,
                                    const stk::mesh::Entity& element, const MDArray& parametric_coordinates, double time_value_optional) 
    {
      PRINT("tmp srk StringFunction::operator(element) getName()= " << getName() << " input_phy_points= " << input_phy_points << " output_values= " << output_values);

      argsAreValid(input_phy_points, output_values);
      argsAreValid(parametric_coordinates, output_values);

      m_element = &element;
      m_have_element = true;
      VERIFY_OP(parametric_coordinates.rank(), ==, 2, "StringFunction::operator() parametric_coordinates rank bad");
      m_parametric_coordinates = MDArray(1, parametric_coordinates.dimension(1));
      int nPoints=parametric_coordinates.dimension(0);
      int spaceDim = parametric_coordinates.dimension(1);
      for (int iPoint = 0; iPoint < nPoints; iPoint++)
        {
          for (int iSpace=0; iSpace < spaceDim; iSpace++)
            m_parametric_coordinates(0, iSpace)=parametric_coordinates(iPoint, iSpace);
          (*this)(input_phy_points, output_values, time_value_optional);
        }

      // reset this else we won't be able to reuse this object correctly
      m_have_element = false;
      m_element = 0;

      // *  Dimensions of parametric_coordinates are required to be ([P],[D])
      //FieldFunction:: void operator()(const stk::mesh::Entity *element, const MDArray& parametric_coordinates, MDArray& out);

    }

    void FieldFunction::operator()(MDArray& input_phy_points, MDArray& output_field_values,
                                   const stk_classic::mesh::Bucket& bucket, const MDArray& parametric_coordinates, double time_value_optional)
    {
      EXCEPTWATCH;
#ifndef NDEBUG
      int num_elements_in_bucket   = bucket.size();
      VERIFY_OP(input_phy_points.dimension(0), ==, num_elements_in_bucket, "FieldFunction::operator() mismatch in input_phy_points and num_elements_in_bucket");
      VERIFY_OP(output_field_values.dimension(0), ==, num_elements_in_bucket, "FieldFunction::operator() mismatch in input_phy_points and num_elements_in_bucket");
#endif
      helper(input_phy_points, output_field_values, bucket, parametric_coordinates, time_value_optional);
    }

 static inline void first_dimensions(MDArray& arr, int arr_offset, int *n_points, int max_rank=3)
 {
   for (int ii = 0; ii < max_rank; ii++)
     {
       n_points[ii] = 1;
     }
   for (int ii = 0; ii < arr_offset; ii++)
     {
       n_points[ii] = arr.dimension(ii);
     }
 }

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 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) );
}