C++ outer函数代码示例

FiniteVolumeEquation<Vector2D> laplacian(Scalar gamma, VectorFiniteVolumeField &phi, Scalar theta)
{
    FiniteVolumeEquation<Vector2D> eqn(phi);
    const VectorFiniteVolumeField &phi0 = phi.oldField(0);

    for (const Cell &cell: phi.cells())
    {
        for (const InteriorLink &nb: cell.neighbours())
        {
            Scalar coeff = gamma * dot(nb.rCellVec(), nb.outwardNorm()) / nb.rCellVec().magSqr();
            eqn.add(cell, nb.cell(), theta * coeff);
            eqn.add(cell, cell, theta * -coeff);
            eqn.addSource(cell, (1. - theta) * coeff * (phi0(nb.cell()) - phi0(cell)));
        }

        for (const BoundaryLink &bd: cell.boundaries())
        {
            Scalar coeff = gamma * dot(bd.rFaceVec(), bd.outwardNorm()) / bd.rFaceVec().magSqr();

            switch (phi.boundaryType(bd.face()))
            {
            case VectorFiniteVolumeField::FIXED:
                eqn.add(cell, cell, theta * -coeff);
                eqn.addSource(cell, theta * coeff * phi(bd.face()));
                eqn.addSource(cell, (1. - theta) * coeff * (phi0(bd.face()) - phi0(cell)));
                break;

            case VectorFiniteVolumeField::NORMAL_GRADIENT:
                break;
            case VectorFiniteVolumeField::SYMMETRY:
            {
                Vector2D tw = bd.outwardNorm().tangentVec().unitVec();

                eqn.add(cell, cell, theta * -coeff);
                eqn.add(cell, cell, theta * coeff * outer(tw, tw));
                eqn.addSource(cell, (1. - theta) * coeff * (dot(phi0(cell), tw) * tw - phi0(cell)));
            }
                break;

            case VectorFiniteVolumeField::PARTIAL_SLIP:
            {
                Vector2D tw = bd.outwardNorm().tangentVec().unitVec();
                Scalar lambda = phi.boundaryRefValue(bd.face()).x;

                Scalar a = lambda != 0. ? lambda * coeff / (lambda * coeff - 1.) : 0.;

                eqn.add(cell, cell, theta * -coeff);
                eqn.add(cell, cell, theta * a * coeff * outer(tw, tw));
                eqn.addSource(cell, (1. - theta) * coeff * (a * dot(phi0(cell), tw) * tw - phi0(cell)));
            }

            default:
                throw Exception("fv", "laplacian<Vector2D>", "unrecognized or unspecified boundary type.");
            }
        }
    }

    return eqn;
}

void CalcDensity(particle_t *SPH){
	#pragma omp parallel for
	for(int i = 0 ; i < N_SPHP ; ++ i){
		SPH[i].div_v = 0;
		SPH[i].omega = 0;
		SPH[i].rot_v = vec3<double>(0, 0, 0);
		SPH[i].p_smth = 0;
		for(int k = 0 ; k < SPH[i].n_ngb ; ++ k){
			int j = SPH[i].ngb_hash[k];
			SPH[i].p_smth += SPH[j].Y * kernel.W(SPH[i].r - SPH[j].r, SPH[i].h);
		}
		//Pressure
		SPH[i].p   = Pressure  (SPH[i].m * SPH[i].p_smth / SPH[i].Y, SPH[i].u);
		SPH[i].c   = SoundSpeed(SPH[i].m * SPH[i].p_smth / SPH[i].Y, SPH[i].u);
		SPH[i].rho = SPH[i].m * SPH[i].p_smth / SPH[i].Y;
		//div v and rot v
		for(int k = 0 ; k < SPH[i].n_ngb ; ++ k){
			int j = SPH[i].ngb_hash[k];
			if(i == j) continue;
			vec3<double> dr = SPH[i].r - SPH[j].r;
			vec3<double> dv = SPH[i].v - SPH[j].v;
			vec3<double> grad_kernel = kernel.gradW(dr, SPH[i].h);
			SPH[i].omega += - SPH[j].Y * (N_DIM / SPH[i].h * kernel.W(dr, SPH[i].h) - abs(dr) / SPH[i].h * abs(grad_kernel));
			double volume = SPH[j].Y / SPH[i].p_smth;
			SPH[i].div_v += - volume * inner(dv, grad_kernel);
			SPH[i].rot_v += - volume * outer(dv, grad_kernel);
		}
		//DEBUG
		SPH[i].er = (SPH[i].p_smth + 7.15 * 3309.0) / (SPH[i].rho * 6.15) / SPH[i].u - 1.0;
		//Balsara 1989
		SPH[i].f = abs(SPH[i].div_v) / (abs(SPH[i].div_v) + abs(SPH[i].rot_v) + 1.0e-4 * SPH[i].c / SPH[i].h);
		//Hosono+ (?)
		SPH[i].omega = 1.0 + SPH[i].h / (N_DIM * SPH[i].p_smth) * SPH[i].omega;
	}
}

void CompoundBody::ComputeInertia()
{
	//http://en.wikipedia.org/wiki/Parallel_axis_theorem
	//Uses the parallel axis theorem to compute a combined
	//Inertia tensor in world space around the center of mass
	//for the entire compound body

	//Skew-symmetric matrices are used as a supplement to
	//cross products.

	AglMatrix inertia = AglMatrix::zeroMatrix();
	AglMatrix id = AglMatrix::identityMatrix();
	for (unsigned int i = 0; i < mChildren.size(); i++)
	{
		AglMatrix orientation = mChildren[i]->GetWorld();
		orientation[15] = 1;
		orientation[14] = 0;
		orientation[13] = 0;
		orientation[12] = 0;
		AglMatrix InertiaWorld = AglMatrix::transpose(orientation) * mChildren[i]->GetLocalInertia() * orientation;
		AglVector3 r = mChildren[i]->GetLocalCenterOfMass() - GetCenterOfMass();

		//OuterProduct
		AglMatrix outer(r.x * r.x, r.x * r.y, r.x * r.z, 0, r.y * r.x, r.y * r.y, r.z * r.x, 0,
						r.z * r.x, r.z * r.y, r.z * r.x, 0, 0, 0, 0, 0);

		inertia += InertiaWorld + (id*AglVector3::dotProduct(r, r) - outer) * mChildren[i]->GetLocalMass();
	}
	mInertiaWorld = inertia;
	mInverseInertiaWorld = AglMatrix::inverse(mInertiaWorld);
}

void TriangulationCDTWindow::OneNestedPolygon()
{
    int numPoints = 7;
    mPoints.resize(numPoints);
    mPoints[0] = { 128.0f, 256.0f };
    mPoints[1] = { 256.0f, 128.0f };
    mPoints[2] = { 384.0f, 256.0f };
    mPoints[3] = { 256.0f, 384.0f };
    mPoints[4] = { 320.0f, 256.0f };
    mPoints[5] = { 256.0f, 192.0f };
    mPoints[6] = { 256.0f, 320.0f };

    std::vector<int> outer(4);
    outer[0] = 0;
    outer[1] = 1;
    outer[2] = 2;
    outer[3] = 3;

    std::vector<int> inner(3);
    inner[0] = 4;
    inner[1] = 5;
    inner[2] = 6;

    Triangulator::Polygon outerPoly = { (int)outer.size(), &outer[0] };
    Triangulator::Polygon innerPoly = { (int)inner.size(), &inner[0] };

    mTriangulator = std::make_unique<Triangulator>(numPoints, &mPoints[0]);
    (*mTriangulator)(outerPoly, innerPoly);
    auto const& triangles = mTriangulator->GetAllTriangles();
    int numTriangles = (int)triangles.size();
    mPMesher = std::make_unique<PlanarMesher>(numPoints, &mPoints[0],
        numTriangles, (int const*)&triangles[0]);

    DrawTriangulation();
}

BaseIF* makeVane(const Real&     thick,
                 const RealVect& normal,
                 const Real&     innerRadius,
                 const Real&     outerRadius,
                 const Real&     offset,
                 const Real&     height)
{
  RealVect zero(D_DECL(0.0,0.0,0.0));
  RealVect xAxis(D_DECL(1.0,0.0,0.0));
  bool inside = true;

  Vector<BaseIF*> vaneParts;

  // Each side of the vane (infinite)
  RealVect normal1(normal);
  RealVect point1(D_DECL(offset+height/2.0,-thick/2.0,0.0));
  PlaneIF plane1(normal1,point1,inside);

  vaneParts.push_back(&plane1);

  RealVect normal2(-normal);
  RealVect point2(D_DECL(offset+height/2.0,thick/2.0,0.0));
  PlaneIF plane2(normal2,point2,inside);

  vaneParts.push_back(&plane2);

  // Make sure we only get something to the right of the origin
  RealVect normal3(D_DECL(0.0,0.0,1.0));
  RealVect point3(D_DECL(0.0,0.0,0.0));
  PlaneIF plane3(normal3,point3,inside);

  vaneParts.push_back(&plane3);

  // Cut off the top and bottom
  RealVect normal4(D_DECL(1.0,0.0,0.0));
  RealVect point4(D_DECL(offset,0.0,0.0));
  PlaneIF plane4(normal4,point4,inside);

  vaneParts.push_back(&plane4);

  RealVect normal5(D_DECL(-1.0,0.0,0.0));
  RealVect point5(D_DECL(offset+height,0.0,0.0));
  PlaneIF plane5(normal5,point5,inside);

  vaneParts.push_back(&plane5);

  // The outside of the inner cylinder
  TiltedCylinderIF inner(innerRadius,xAxis,zero,!inside);

  vaneParts.push_back(&inner);

  // The inside of the outer cylinder
  TiltedCylinderIF outer(outerRadius,xAxis,zero,inside);

  vaneParts.push_back(&outer);

  IntersectionIF* vane = new IntersectionIF(vaneParts);

  return vane;
}

int main()
{
    int ret;

    // OK
    ret = foo();
    if (ret < 0)
    {
        xmlSecInternalError("foo", NULL);
    }

    // OK
    ret = outer(1, strlen("x"));
    if (ret < 0)
    {
        xmlSecInternalError("outer", NULL);
    }

    // KO
    ret = foo();
    if (ret < 0)
    {
        xmlSecInternalError("bar", NULL);
    }
}

	sf::ConvexShape star(unsigned int nbStarPoints, float innerRadius, float outerRadius,
		const sf::Color& fillColor, float outlineThickness, const sf::Color& outlineColor)
	{
		assert(innerRadius > 0.f);
		assert(outerRadius > innerRadius);
		assert(outlineThickness >= 0.f);

		// Calculate points of the inner, regular polygon and the outer star points
		PolarVector2f      inner(innerRadius, 0.f);
		PolarVector2f      outer(outerRadius, 0.f);
	
		sf::ConvexShape shape;
		shape.setFillColor(fillColor);
		shape.setOutlineThickness(outlineThickness);
		shape.setOutlineColor(outlineColor);
		
		// Step around and alternately add inner and outer points
		for (unsigned int points = 0; points < nbStarPoints; ++points)
		{
			inner.phi = 360.f * points / nbStarPoints;
			outer.phi = inner.phi + 180.f / nbStarPoints;
		
			addPoint(shape, inner);
			addPoint(shape, outer);
		}

		return shape;
	}

sk_sp<SkSpecialImage> SkComposeImageFilter::onFilterImage(SkSpecialImage* source,
                                                          const Context& ctx,
                                                          SkIPoint* offset) const {
    // The bounds passed to the inner filter must be filtered by the outer
    // filter, so that the inner filter produces the pixels that the outer
    // filter requires as input. This matters if the outer filter moves pixels.
    SkIRect innerClipBounds;
    innerClipBounds = this->getInput(0)->filterBounds(ctx.clipBounds(), ctx.ctm());
    Context innerContext(ctx.ctm(), innerClipBounds, ctx.cache());
    SkIPoint innerOffset = SkIPoint::Make(0, 0);
    sk_sp<SkSpecialImage> inner(this->filterInput(1, source, innerContext, &innerOffset));
    if (!inner) {
        return nullptr;
    }

    SkMatrix outerMatrix(ctx.ctm());
    outerMatrix.postTranslate(SkIntToScalar(-innerOffset.x()), SkIntToScalar(-innerOffset.y()));
    SkIRect clipBounds = ctx.clipBounds();
    clipBounds.offset(-innerOffset.x(), -innerOffset.y());
    Context outerContext(outerMatrix, clipBounds, ctx.cache());

    SkIPoint outerOffset = SkIPoint::Make(0, 0);
    sk_sp<SkSpecialImage> outer(this->filterInput(0, inner.get(), outerContext, &outerOffset));
    if (!outer) {
        return nullptr;
    }

    *offset = innerOffset + outerOffset;
    return outer;
}

//---------------------------------------------------------
void xyztorst
(
  const DVec& X,  // [in]
  const DVec& Y,  // [in]
  const DVec& Z,  // [in]
        DVec& r,  // [out]
        DVec& s,  // [out]
        DVec& t   // [out]
)
//---------------------------------------------------------
{
  // function [r,s,t] = xyztorst(x, y, z)
  // Purpose : Transfer from (x,y,z) in equilateral tetrahedron
  //           to (r,s,t) coordinates in standard tetrahedron

  double sqrt3=sqrt(3.0), sqrt6=sqrt(6.0); int Nc=X.size();
  DVec v1(3),v2(3),v3(3),v4(3);
  DMat tmat1(3,Nc), A(3,3), rhs; 
  
  v1(1)=(-1.0);  v1(2)=(-1.0/sqrt3);  v1(3)=(-1.0/sqrt6);
  v2(1)=( 1.0);  v2(2)=(-1.0/sqrt3);  v2(3)=(-1.0/sqrt6);
  v3(1)=( 0.0);  v3(2)=( 2.0/sqrt3);  v3(3)=(-1.0/sqrt6);
  v4(1)=( 0.0);  v4(2)=( 0.0      );  v4(3)=( 3.0/sqrt6);

  // back out right tet nodes
  tmat1.set_row(1,X); tmat1.set_row(2,Y); tmat1.set_row(3,Z);
  rhs = tmat1 - 0.5*outer(v2+v3+v4-v1, ones(Nc));
  A.set_col(1,0.5*(v2-v1)); A.set_col(2,0.5*(v3-v1)); A.set_col(3,0.5*(v4-v1));

  DMat RST = A|rhs;

  r=RST.get_row(1); s=RST.get_row(2); t=RST.get_row(3);
}

 void update_derivs(const vector<int>& toCoords)
 {
     const vector<int>* fromCoords = add_delay(toCoords);
     if (fromCoords) 
     {
         outer(this->from->out_acts(*fromCoords), derivs().begin(), this->to->inputErrors[toCoords]);
     }
 }

gmMatrix3 CSGUnion::hess(const gmVector3 & x)
{
  if ((m_f!=NULL) && (m_g!=NULL))
    {
      double fx = m_f->proc(x); 
      double gx = m_g->proc(x);
      gmVector3 dfx = m_f->grad(x);
      gmVector3 dgx = m_g->grad(x);

      return hff(fx,gx) * outer(dfx,dfx) + 
             hfg(fx,gx) * outer(dfx,dgx) +
             hfg(fx,gx) * outer(dgx,dfx) + 
             hgg(fx,gx) * outer(dgx,dgx) +
             hf(fx,gx)  * m_f->hess(x)   + 
             hg(fx,gx)  * m_g->hess(x);
    }
  else
    return gmMatrix3();
}

void CalcConservativeVaruables(const particle_t *SPH, system_t *system){
	system->energy = 0;
	system->linear_momentum  = vec3<double>(0, 0, 0);
	system->angular_momentum = vec3<double>(0, 0, 0);
	#pragma omp parallel for
	for(int i = 0 ; i < N_SPHP ; ++ i){
		system->linear_momentum  += SPH[i].m * SPH[i].v;
		system->angular_momentum += SPH[i].m * outer(SPH[i].r, SPH[i].v);
		system->energy           += SPH[i].m * (0.5 * abs2(SPH[i].v) + SPH[i].u);
	}
}

int main (void) {
  char *dummy = __builtin_alloca(alloca_size());
  int retval;

  fprintf(stderr, "main: dummy = %p\n", dummy);
  retval = outer(dummy);
  fprintf(stderr, "main: outer returned %d\n", retval);
  if (retval == 0)
    abort();
  return 0;
}

  namespace numpy
  {
    template <class T0, size_t N0, class T1, size_t N1>
    types::ndarray<decltype(std::declval<T0>() + std::declval<T1>()), 2>
    outer(types::ndarray<T0, N0> const &a, types::ndarray<T1, N1> const &b);

    template <class T0, size_t N0, class E1>
    auto outer(types::ndarray<T0, N0> const &a, E1 const &b)
        -> decltype(outer(a, asarray(b)));

    template <class E0, class T1, size_t N1>
    auto outer(E0 const &a, types::ndarray<T1, N1> const &b)
        -> decltype(outer(asarray(a), b));

    template <class E0, class E1>
    auto outer(E0 const &a, E1 const &b)
        -> decltype(outer(asarray(a), asarray(b)));

    PROXY_DECL(pythonic::numpy, outer);
  }

const char *ScriptObject::GetFullName()
{
	if( object_class() && outer() )
	{
		static std::string full_name;
		full_name = GetName();

		for( ScriptObject *object_outer = outer(); object_outer != NULL; object_outer = object_outer->outer() )
		{
			full_name.insert( 0, "." );
			full_name.insert( 0, object_outer->GetName() );
		}

		full_name.insert( 0, " " );
		full_name.insert( 0, object_class()->GetName() );

		return full_name.c_str();
	}

	return "Failed to get name";
}