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

int main(int argc, char* argv[])
{
  srand((unsigned)time(0));
  Vector dVec(5);
  Matrix m(3, 3);
  m(0, 0) = -1.0; m(0, 1) = 1.0; m(0, 2) = 2.0;
  m(1, 0) = 2.0; m(1, 1) = 3.0; m(1, 2) = 3.0;
  m(2, 0) = 4.0; m(2, 1) = 2.0; m(2, 2) = 1.0;
  m.print();
  dVec.resize(3);
  dVec[0] = 0.78045432; dVec[1] = 0.27960367; dVec[2] = 0.55920734;
  Vector d2(3);
  d2 = dVec*m;
  d2.print();
  std::cout << "\n\n\n";
  d2[0] = 3.5765; d2[1] = 2.7377; d2[2] = 2.9590;
  m = outer(dVec, d2);
  std::cout << "\n\n outer product \n\n";
  m.print();
  std::cout << "\n\n\n";
  dVec.print();

  dVec.resize(20);
  for (int i = 0; i < 20; i++){
    dVec[i] = rand()%100 + 1;
  }
  dVec.print();
  dVec.sort();
  dVec.print();
  int dim = 5;
  double PRECISION = 1e-12;
  Vector v(dim);
  v[0] = 0.0; v[1] = 5.4; v[2] = 3.2; v[3] = 0.0; v[4] = 4.1;
  Matrix P(dim, dim, 0.0);
  for (int j = 0; j < dim; j++) { P(j, j) = 1.0; } // Turn into identity                                                                                         
  int k= dim-1; // Keep track of last non-zero entry                                                                                                            
  for (int j = dim-1; j > -1; j--){ // Start at bottom and work up                                                                   
    if (fabs(v(j))< PRECISION && k!=j ){
      std::cout << "\n Found a zero. \n";
      // And swap the two rows                                                                                                                                   
      P.swapRows(j, k);
      k--;
    }
  }
  std::cout << "\n\n";
  P.print();
  v = P*v;
  std::cout << "\n\n";
  v.print();

  return 0;
}

int main(int argc, char *argv[])
{

  int i = 0;
  char *ret;
  for( ;i < 10; i++){
    printf("%d: %s\n", i, ret);
    ret = outer();
  }
	

  return 0;
}

// Helper function for managing labels and their target addresses.
// Returns a sensible address, and if it is not the label's final
// address, notes the dependency (at 'branch_pc') on the label.
address CodeSection::target(Label& L, address branch_pc) {
  if (L.is_bound()) {
    int loc = L.loc();
    if (index() == CodeBuffer::locator_sect(loc)) {
      return start() + CodeBuffer::locator_pos(loc);
    } else {
      return outer()->locator_address(loc);
    }
  } else {
    assert(allocates2(branch_pc), "sanity");
    address base = start();
    int patch_loc = CodeBuffer::locator(branch_pc - base, index());
    L.add_patch_at(outer(), patch_loc);

    // Need to return a pc, doesn't matter what it is since it will be
    // replaced during resolution later.
    // Don't return NULL or badAddress, since branches shouldn't overflow.
    // Don't return base either because that could overflow displacements
    // for shorter branches.  It will get checked when bound.
    return branch_pc;
  }
}

int main ()
{
	std::list<int> t, l{4, 3, 2, 5, 0, 1};
	std::cout << l << std::endl;

	auto b = l.begin(), e = l.begin();
	std::advance(b, 1);
	std::advance(e, 4);
	bool in = true;

	in ? inner(b, e) : outer(l, b, e);
	std::cout << l << std::endl;
}

NumericMatrix polyOuter(const NumericMatrix &Thetas, const vector<double> &Pk,
	const vector<double> &Pk_1, const vector<double> &PQ_1, const vector<double> &PQ,
	const vector<double> &dif1sq, const vector<double> &dif1)
{
	const int nfact = Thetas.ncol();
	NumericMatrix d2Louter(nfact,nfact), outer(nfact,nfact);
	vector<double> temp(nfact);

	for(int n = 0; n < Thetas.nrow(); ++n){
		for(int i = 0; i < nfact; ++i)
			for(int j = 0; j < nfact; ++j)
				outer(i,j) = Thetas(n,i) * Thetas(n,j);
		for(int i = 0; i < nfact; ++i)
			temp[i] =  (PQ_1[n] * Thetas(n,i) - PQ[n] * Thetas(n,i));
		for(int i = 0; i < nfact; ++i)
			for(int j = 0; j < nfact; ++j)
				d2Louter(i,j) += (-1) * dif1sq[n] * temp[i] * temp[j] +
				    (dif1[n] * (Pk_1[n] * (1.0 - Pk_1[n]) * (1.0 - 2.0 * Pk_1[n]) *
				    outer(i,j) - Pk[n] * (1.0 - Pk[n]) * (1.0 - 2.0 * Pk[n]) * outer(i,j)));
	}
	return d2Louter;
}