C++ GridFunction类代码示例

arma::Mat<ResultType>
ElementaryPotentialOperator<BasisFunctionType, KernelType, ResultType>::
    evaluateAtPoints(
        const GridFunction<BasisFunctionType, ResultType> &argument,
        const arma::Mat<CoordinateType> &evaluationPoints,
        const QuadratureStrategy &quadStrategy,
        const EvaluationOptions &options) const {
  if (evaluationPoints.n_rows != argument.grid()->dimWorld())
    throw std::invalid_argument(
        "ElementaryPotentialOperator::evaluateAtPoints(): "
        "the number of coordinates of each evaluation point must be "
        "equal to the dimension of the space containing the surface "
        "on which the grid function 'argument' is defined");

  if (options.evaluationMode() == EvaluationOptions::DENSE) {
    std::unique_ptr<Evaluator> evaluator =
        makeEvaluator(argument, quadStrategy, options);

    // right now we don't bother about far and near field
    // (this might depend on evaluation options)
    arma::Mat<ResultType> result;
    evaluator->evaluate(Evaluator::FAR_FIELD, evaluationPoints, result);
    return result;
  } else if (options.evaluationMode() == EvaluationOptions::ACA) {
    AssembledPotentialOperator<BasisFunctionType, ResultType> assembledOp =
        assemble(argument.space(), make_shared_from_ref(evaluationPoints),
                 quadStrategy, options);
    return assembledOp.apply(argument);
  } else
    throw std::invalid_argument(
        "ElementaryPotentialOperator::evaluateAtPoints(): "
        "Invalid evaluation mode");
}

void ParGridFunction::ComputeFlux(
   BilinearFormIntegrator &blfi,
   GridFunction &flux, int wcoef, int subdomain)
{
   ParFiniteElementSpace *ffes =
      dynamic_cast<ParFiniteElementSpace*>(flux.FESpace());
   MFEM_VERIFY(ffes, "the flux FE space must be ParFiniteElementSpace");

   Array<int> count(flux.Size());
   SumFluxAndCount(blfi, flux, count, wcoef, subdomain);

   // Accumulate flux and counts in parallel
   ffes->GroupComm().Reduce<double>(flux, GroupCommunicator::Sum);
   ffes->GroupComm().Bcast<double>(flux);

   ffes->GroupComm().Reduce<int>(count, GroupCommunicator::Sum);
   ffes->GroupComm().Bcast<int>(count);

   // complete averaging
   for (int i = 0; i < count.Size(); i++)
   {
      if (count[i] != 0) { flux(i) /= count[i]; }
   }

   if (ffes->Nonconforming())
   {
      // On a partially conforming flux space, project on the conforming space.
      // Using this code may lead to worse refinements in ex6, so we do not use
      // it by default.

      // Vector conf_flux;
      // flux.ConformingProject(conf_flux);
      // flux.ConformingProlongate(conf_flux);
   }
}

VisualizationSceneVector3d::VisualizationSceneVector3d(GridFunction &vgf)
{
   FiniteElementSpace *fes = vgf.FESpace();
   if (fes == NULL || fes->GetVDim() != 3)
   {
      cout << "VisualizationSceneVector3d::VisualizationSceneVector3d" << endl;
      exit(1);
   }

   VecGridF = &vgf;

   mesh = fes->GetMesh();

   sfes = new FiniteElementSpace(mesh, fes->FEColl(), 1, fes->GetOrdering());
   GridF = new GridFunction(sfes);

   solx = new Vector(mesh->GetNV());
   soly = new Vector(mesh->GetNV());
   solz = new Vector(mesh->GetNV());

   vgf.GetNodalValues(*solx, 1);
   vgf.GetNodalValues(*soly, 2);
   vgf.GetNodalValues(*solz, 3);

   sol = new Vector(mesh->GetNV());

   Init();
}

void Computation::computeNewVelocities(GridFunction& u, GridFunction& v,
		GridFunction& f, GridFunction& g, GridFunction& p, RealType deltaT) {
	GridFunction branch_1(p.griddimension);
	MultiIndexType begin, end;
	begin[0] = 1;
	end[0] = u.griddimension[0] - 3;
	begin[1] = 1;
	end[1] = u.griddimension[1] - 2;
	// u Update
	PointType delta;
	delta[0] = SimIO.para.deltaX;
	delta[1] = SimIO.para.deltaY;
	Px(branch_1, p, delta);
	f.AddToGridFunction(begin, end, -1.0 * deltaT, branch_1);
	u.SetGridFunction(begin, end, 1.0, f);
	// v Update
	begin[0] = 1;
	end[0] = v.griddimension[0] - 2;
	begin[1] = 1;
	end[1] = v.griddimension[1] - 3;
	GridFunction branch_2(p.griddimension);
	Py(branch_2, p, delta);
	g.AddToGridFunction(begin, end, -1.0 * deltaT, branch_2);
	v.SetGridFunction(begin, end, 1.0, g);

}

void VVy(GridFunction& output, GridFunction& v, const RealType alpha,
		const PointType& h) {

	MultiIndexType begin, end;

	begin[0] = 1;
	end[0] = v.griddimension[0] - 2;
	begin[1] = 1;
	end[1] = v.griddimension[1] - 3;

	Stencil stencil_1(3, h);
	stencil_1.setVVy_1Stencil();
	stencil_1.ApplyStencilOperator(begin, end, begin, end, v, output);
	output.MultiplyGridFunctions(begin, end, output);

	GridFunction branch_2(v.griddimension);
	Stencil stencil_2(3, h);
	stencil_2.setVVy_2Stencil();
	stencil_2.ApplyStencilOperator(begin, end, begin, end, v, branch_2);
	branch_2.MultiplyGridFunctions(begin, end, branch_2);

	output.AddToGridFunction(begin, end, -1.0, branch_2);
	output.ScaleGridFunction(begin, end, 1.0 / h[1]);

	GridFunction branch_3(v.griddimension);
	Stencil stencil_3(3, h);
	stencil_3.setVVy_3Stencil();
	stencil_3.ApplyStencilOperator(begin, end, begin, end, v, branch_3);
	branch_3.MultiplyGridFunctions(begin, end, branch_3);

	GridFunction branch_4(v.griddimension);
	Stencil stencil_4(3, h);
	stencil_4.setVVy_4Stencil();
	stencil_4.ApplyStencilOperator(begin, end, begin, end, v, branch_4);
	branch_4.MultiplyGridFunctions(begin, end, branch_4);

	branch_3.MultiplyGridFunctions(begin, end, branch_4);

	GridFunction branch_5(v.griddimension);
	Stencil stencil_5(3, h);
	stencil_5.setVVy_5Stencil();
	stencil_5.ApplyStencilOperator(begin, end, begin, end, v, branch_5);
	branch_5.MultiplyGridFunctions(begin, end, branch_5);

	GridFunction branch_6(v.griddimension);
	Stencil stencil_6(3, h);
	stencil_6.setVVy_6Stencil();
	stencil_6.ApplyStencilOperator(begin, end, begin, end, v, branch_6);
	branch_6.MultiplyGridFunctions(begin, end, branch_6);

	branch_5.MultiplyGridFunctions(begin, end, branch_6);

	branch_3.AddToGridFunction(begin, end, -1.0, branch_5);

	branch_3.ScaleGridFunction(begin, end, alpha / h[1]);

	output.AddToGridFunction(begin, end, 1.0, branch_3);

}

void UVx(GridFunction& output, GridFunction& u, GridFunction& v,
		const RealType alpha, const PointType& h) {
	MultiIndexType begin, end;

	begin[0] = 1;
	end[0] = u.griddimension[0] - 2;
	begin[1] = 1;
	end[1] = u.griddimension[1] - 3;

	Stencil stencil_1(3, h);
	stencil_1.setUVx_1Stencil();
	stencil_1.ApplyStencilOperator(begin, end, begin, end, u, output);

	GridFunction branch_2(u.griddimension);
	Stencil stencil_2(3, h);
	stencil_2.setUVx_2Stencil();
	stencil_2.ApplyStencilOperator(begin, end, begin, end, v, branch_2);
	output.MultiplyGridFunctions(begin, end, branch_2);

	GridFunction branch_3(u.griddimension);
	Stencil stencil_3(3, h);
	stencil_3.setUVx_3Stencil();
	stencil_3.ApplyStencilOperator(begin, end, begin, end, u, branch_3);

	GridFunction branch_4(u.griddimension);
	Stencil stencil_4(3, h);
	stencil_4.setUVx_4Stencil();
	stencil_4.ApplyStencilOperator(begin, end, begin, end, v, branch_4);
	branch_3.MultiplyGridFunctions(begin, end, branch_4);
	output.AddToGridFunction(begin, end, -1.0, branch_3);
	output.ScaleGridFunction(begin, end, 1.0 / h[0]);
//bis hier 1. Zeile
	GridFunction branch_5(u.griddimension);
	Stencil stencil_5(3, h);
	stencil_5.setUVx_5Stencil();
	stencil_5.ApplyStencilOperator(begin, end, begin, end, u, branch_5);

	GridFunction branch_6(u.griddimension);
	Stencil stencil_6(3, h);
	stencil_6.setUVx_6Stencil();
	stencil_6.ApplyStencilOperator(begin, end, begin, end, v, branch_6);
	branch_5.MultiplyGridFunctions(begin, end, branch_6);

	GridFunction branch_7(u.griddimension);
	Stencil stencil_7(3, h);
	stencil_7.setUVx_7Stencil();
	stencil_7.ApplyStencilOperator(begin, end, begin, end, u, branch_7);

	GridFunction branch_8(u.griddimension);
	Stencil stencil_8(3, h);
	stencil_8.setUVx_8Stencil();
	stencil_8.ApplyStencilOperator(begin, end, begin, end, v, branch_8);
	branch_7.MultiplyGridFunctions(begin, end, branch_8);
	branch_5.AddToGridFunction(begin, end, -1.0, branch_7);
	branch_5.ScaleGridFunction(begin, end, alpha * 1.0 / h[0]);
//bis hier 2. Zeile
	output.AddToGridFunction(begin, end, 1.0, branch_5);
}

void Computation::ComputeHeatfunction(GridFunction& h, GridFunction& t, GridFunction& u, RealType deltaT) {
	for (int i = 0; i <= h.griddimension[0]-2; i++){
		for (int j = 1; j <= h.griddimension[1]-2; j++){
			h.getGridFunction()[i][j] = h.getGridFunction()[i][j-1]
			                + deltaT * (SimIO.para.re * SimIO.para.Pr * u.getGridFunction()[i][j]
			                  * (t.getGridFunction()[i+1][j] + t.getGridFunction()[i][j])/2.0
			                  - (t.getGridFunction()[i+1][j] - t.getGridFunction()[i][j])/ SimIO.para.deltaX);
		}
	}

}

std::auto_ptr<typename ElementaryPotentialOperator<
BasisFunctionType, KernelType, ResultType>::Evaluator>
ElementaryPotentialOperator<BasisFunctionType, KernelType, ResultType>::
makeEvaluator(
        const GridFunction<BasisFunctionType, ResultType>& argument,
        const QuadratureStrategy& quadStrategy,
        const EvaluationOptions& options) const
{
    // Collect the standard set of data necessary for construction of
    // evaluators and assemblers
    typedef Fiber::RawGridGeometry<CoordinateType> RawGridGeometry;
    typedef std::vector<const Fiber::Shapeset<BasisFunctionType>*> ShapesetPtrVector;
    typedef std::vector<std::vector<ResultType> > CoefficientsVector;
    typedef LocalAssemblerConstructionHelper Helper;

    shared_ptr<RawGridGeometry> rawGeometry;
    shared_ptr<GeometryFactory> geometryFactory;
    shared_ptr<Fiber::OpenClHandler> openClHandler;
    shared_ptr<ShapesetPtrVector> shapesets;

    const Space<BasisFunctionType>& space = *argument.space();
    shared_ptr<const Grid> grid = space.grid();
    Helper::collectGridData(space,
                            rawGeometry, geometryFactory);
    Helper::makeOpenClHandler(options.parallelizationOptions().openClOptions(),
                              rawGeometry, openClHandler);
    Helper::collectShapesets(space, shapesets);

    // In addition, get coefficients of argument's expansion in each element
    const GridView& view = space.gridView();
    const int elementCount = view.entityCount(0);

    shared_ptr<CoefficientsVector> localCoefficients =
            boost::make_shared<CoefficientsVector>(elementCount);

    std::auto_ptr<EntityIterator<0> > it = view.entityIterator<0>();
    for (int i = 0; i < elementCount; ++i) {
        const Entity<0>& element = it->entity();
        argument.getLocalCoefficients(element, (*localCoefficients)[i]);
        it->next();
    }

    // Now create the evaluator
    return quadStrategy.makeEvaluatorForIntegralOperators(
                geometryFactory, rawGeometry,
                shapesets,
                make_shared_from_ref(kernels()),
                make_shared_from_ref(trialTransformations()),
                make_shared_from_ref(integral()),
                localCoefficients,
                openClHandler,
                options.parallelizationOptions());
}

void Computation::setBoundaryG(GridFunction& g, GridFunction& v) {

	MultiIndexType begin, end;
	begin[0] = 1;
	end[0] = g.griddimension[0] - 2;
	begin[1] = 0;
	end[1] = 0;
	g.SetGridFunction(begin, end, 1.0, v);

	begin[0] = 1;
	end[0] = g.griddimension[0] - 2;
	begin[1] = g.griddimension[1] - 2;
	end[1] = g.griddimension[1] - 2;
	g.SetGridFunction(begin, end, 1.0, v);
}

void Computation::setBoundaryF(GridFunction& f, GridFunction& u) {
	MultiIndexType begin, end;
	begin[0] = 0;
	end[0] = 0;
	begin[1] = 1;
	end[1] = f.griddimension[1] - 2;
	f.SetGridFunction(begin, end, 1.0, u);

// F_iMax,j=u_iMax+1,j
	begin[0] = f.griddimension[0] - 2;
	end[0] = f.griddimension[0] - 2;
	begin[1] = 1;
	end[1] = f.griddimension[1] - 2;
	f.SetGridFunction(begin, end, 1.0, u);
}

void Stencil::ApplyUVxStencilOperator(const MultiIndexType& gridreadbegin,
		const MultiIndexType& gridreadend,
		const MultiIndexType& gridwritebegin,
		const MultiIndexType& gridwriteend,
		const GridFunctionType& sourcegridfunctionU,
		const GridFunctionType& sourcegridfunctionV,
		GridFunction& imagegridfunction,
		RealType alpha){
	RealType tmp;
	for (IndexType i=gridwritebegin[0]; i<=gridwriteend[0]; i++){
		for (IndexType j=gridwritebegin[1]; j<=gridwriteend[1]; j++){
			tmp = 0.25*((sourcegridfunctionU[i][j]+sourcegridfunctionU[i][j+1])*
						(sourcegridfunctionV[i][j]+sourcegridfunctionV[i+1][j])-
						(sourcegridfunctionU[i-1][j]+sourcegridfunctionU[i-1][j+1])*
						(sourcegridfunctionV[i-1][j]+sourcegridfunctionV[i][j])
			)+
							  alpha*0.25*(abs(sourcegridfunctionU[i][j]+sourcegridfunctionU[i][j+1])*
									  (sourcegridfunctionV[i][j]-sourcegridfunctionV[i+1][j])-
									      abs(sourcegridfunctionU[i-1][j]+sourcegridfunctionU[i-1][j+1])*
									      (sourcegridfunctionV[i-1][j]-sourcegridfunctionV[i][j])
									  );

			imagegridfunction.SetGridFunction(i,j,tmp/h[0]);
		}
	}
}

void Computation::setBoundaryP(GridFunction& p) {
	MultiIndexType begin, end;

	if (SimIO.para.world_rank == 0) {
		// p_0,j = p_1,j
		begin[0] = 0;
		end[0] = 0;
		begin[1] = 1;
		end[1] = p.griddimension[1] - 2;
		MultiIndexType Offset;
		Offset[0] = 1;
		Offset[1] = 0;
		p.SetGridFunction(begin, end, 1.0, Offset);
	}

	if (SimIO.para.world_rank == 1) {
		// p_iMax+1,j = p_iMax,j
		begin[0] = p.griddimension[0] - 1;
		end[0] = p.griddimension[0] - 1;
		begin[1] = 1;
		end[1] = p.griddimension[1] - 2;
		MultiIndexType Offset;
		Offset[0] = -1;
		Offset[1] = 0;
		p.SetGridFunction(begin, end, 1.0, Offset);
	}
	// p_i,0 = p_i,1
	begin[0] = 1;
	end[0] = p.griddimension[0] - 2;
	begin[1] = 0;
	end[1] = 0;
	MultiIndexType Offset;
	Offset[0] = 0;
	Offset[1] = 1;
	p.SetGridFunction(begin, end, 1.0, Offset);

	// p_i,jMax+1 = p_i,jMax
	begin[0] = 1;
	end[0] = p.griddimension[0] - 2;
	begin[1] = p.griddimension[1] - 1;
	end[1] = p.griddimension[1] - 1;
	Offset[0] = 0;
	Offset[1] = -1;
	p.SetGridFunction(begin, end, 1.0, Offset);

}

std::auto_ptr<InterpolatedFunction<ResultType> >
ElementaryPotentialOperator<BasisFunctionType, KernelType, ResultType>::
evaluateOnGrid(
        const GridFunction<BasisFunctionType, ResultType>& argument,
        const Grid& evaluationGrid,
        const QuadratureStrategy& quadStrategy,
        const EvaluationOptions& options) const
{
    if (evaluationGrid.dimWorld() != argument.grid()->dimWorld())
        throw std::invalid_argument(
                "ElementaryPotentialOperator::evaluateOnGrid(): "
                "the evaluation grid and the surface on which the grid "
                "function 'argument' is defined must be embedded in a space "
                "of the same dimension");

    // Get coordinates of interpolation points, i.e. the evaluationGrid's vertices

    std::auto_ptr<GridView> evalView = evaluationGrid.leafView();
    const int evalGridDim = evaluationGrid.dim();
    const int evalPointCount = evalView->entityCount(evalGridDim);
    arma::Mat<CoordinateType> evalPoints(evalGridDim, evalPointCount);

    const IndexSet& evalIndexSet = evalView->indexSet();
    // TODO: extract into template function, perhaps add case evalGridDim == 1
    if (evalGridDim == 2) {
        const int vertexCodim = 2;
        std::auto_ptr<EntityIterator<vertexCodim> > it =
                evalView->entityIterator<vertexCodim>();
        while (!it->finished()) {
            const Entity<vertexCodim>& vertex = it->entity();
            const Geometry& geo = vertex.geometry();
            const int vertexIndex = evalIndexSet.entityIndex(vertex);
            arma::Col<CoordinateType> activeCol(evalPoints.unsafe_col(vertexIndex));
            geo.getCenter(activeCol);
            it->next();
        }
    } else if (evalGridDim == 3) {
        const int vertexCodim = 3;
        std::auto_ptr<EntityIterator<vertexCodim> > it =
                evalView->entityIterator<vertexCodim>();
        while (!it->finished()) {
            const Entity<vertexCodim>& vertex = it->entity();
            const Geometry& geo = vertex.geometry();
            const int vertexIndex = evalIndexSet.entityIndex(vertex);
            arma::Col<CoordinateType> activeCol(evalPoints.unsafe_col(vertexIndex));
            geo.getCenter(activeCol);
            it->next();
        }
    }

    arma::Mat<ResultType> result;
    result = evaluateAtPoints(argument, evalPoints, quadStrategy, options);

    return std::auto_ptr<InterpolatedFunction<ResultType> >(
                new InterpolatedFunction<ResultType>(evaluationGrid, result));
}

void Computation::setBoundaryTN(GridFunction& T, RealType (*TO)(RealType),
		RealType (*TU)(RealType), RealType (*TL)(RealType),
		RealType (*TR)(RealType)) {
	MultiIndexType begin, end;
	// T_i,0

	begin[0] = 1;
	end[0] = T.griddimension[0] - 2;
	begin[1] = 0;
	end[1] = 0;
	MultiIndexType Offset;
	Offset[0] = 0;
	Offset[1] = 1;
	T.SetGridFunction(begin, end, TU, false, SimIO.para.deltaX);
	T.ScaleGridFunction(begin, end, SimIO.para.deltaY);
	T.AddToGridFunction(begin, end, 1.0, T, Offset);

	// T_i,jmax+1

	begin[0] = 1;
	end[0] = T.griddimension[0] - 2;
	begin[1] = T.griddimension[1] - 1;
	end[1] = T.griddimension[1] - 1;
	//MultiIndexType Offset;
	Offset[0] = 0;
	Offset[1] = -1;
	T.SetGridFunction(begin, end, TO, false, SimIO.para.deltaX);
	T.ScaleGridFunction(begin, end, SimIO.para.deltaY);
	T.AddToGridFunction(begin, end, 1.0, T, Offset);
}

void Computation::setBoundaryV(GridFunction& v) {
	MultiIndexType begin, end;

	// v_i,0 = 0
	begin[0] = 1;
	end[0] = v.griddimension[0] - 2;
	begin[1] = 0;
	end[1] = 0;
	v.SetGridFunction(begin, end, 0.0);

	// v_i,jMax =0
	begin[0] = 1;
	end[0] = v.griddimension[0] - 2;
	begin[1] = v.griddimension[1] - 2;
	end[1] = v.griddimension[1] - 2;
	v.SetGridFunction(begin, end, 0.0);
	if (SimIO.para.world_rank == 0) {
		// v_0,j = -v_1,j
		begin[0] = 0;
		end[0] = 0;
		begin[1] = 1;
		end[1] = v.griddimension[1] - 2;
		MultiIndexType Offset;
		Offset[0] = 1;
		Offset[1] = 0;
		v.SetGridFunction(begin, end, -1.0, Offset);
	}
	if (SimIO.para.world_rank == 1) {
		// v_iMax+1,j = -v_iMax,j
		begin[0] = v.griddimension[0] - 1;
		end[0] = v.griddimension[0] - 1;
		begin[1] = 1;
		end[1] = v.griddimension[1] - 2;
		MultiIndexType Offset;
		Offset[0] = -1;
		Offset[1] = 0;
		v.SetGridFunction(begin, end, -1.0, Offset);
	}
}