C++ Face::Bulk方法代码示例

void ConvectionDiffusion::DiffusionFluxes(abstract_variable & param, Face fKL, variable & fluxT, variable & corrK, variable & corrL, variable & corrK_cK, variable & corrL_cL) const
{
	if( !BuildFlux(fKL) ) return;

	get_variable<variable> Func(param);
	//diffusion part
	if( tag_K.isValid() )
	{
		Cell cK = fKL->BackCell();
		Cell cL = fKL->FrontCell();
		//diffusion data
		real       coefT  = fKL.Real           (tag_DIFF_CT);
		real       rhsT   = fKL.Real           (tag_DIFF_RT);
		//corrector in back cell
		real_array coefsK = fKL->RealArray     (tag_DIFF_CB);
		ref_array  elemsK = fKL->ReferenceArray(tag_DIFF_EB);
		real       rhsK   = fKL->Real          (tag_DIFF_RB);
		//corrector in front cell
		real_array coefsL = fKL->RealArray     (tag_DIFF_CF);
		ref_array  elemsL = fKL->ReferenceArray(tag_DIFF_EF);
		real       rhsL   = fKL->Real          (tag_DIFF_RF);
		//flag
		bulk       flag   = fKL->Bulk          (tag_DIFF_F);

		//two-point part
		fluxT -= rhsT;
		if( flag > 1 ) fluxT  -= coefT*(Func(cL) - Func(cK));
		else fluxT  += coefT*Func(cK);

		if( flag > 1 ) //full nonlinear part
		{
			if( perform_correction_diffusion )
			{
				corrK -= rhsK;
				corrL -= rhsL;
				for(enumerator k = 0; k < elemsK.size(); ++k)
					if( elemsK[k].isValid() ) corrK -= Func(elemsK[k])*coefsK[k];
				for(enumerator k = 0; k < elemsL.size(); ++k)
					if( elemsL[k].isValid() ) corrL -= Func(elemsL[k])*coefsL[k];
				//if( scheme_type == NTPFA || scheme_type == NTPFA_PICARD )
				{
					if( elemsK.back().isValid() ) corrK_cK -= Func(elemsK.back())*coefsK.back();
					if( elemsL.back().isValid() ) corrL_cL -= Func(elemsL.back())*coefsL.back();
				}
			}
		}
		else //one-sided corrector
		{
			if( flag == 1 && perform_correction_diffusion ) // only fluxK is present
			{
				corrK -= rhsK;
				for(enumerator k = 0; k < elemsK.size(); ++k)
					if( elemsK[k].isValid() ) corrK -= Func(elemsK[k])*coefsK[k];
			} //no corrections otherwise
		} //flag
	} //diffusion part
}

void ConvectionDiffusion::AdvectionFluxes(abstract_variable & param, Face fKL, variable & fluxT, variable & corrK, variable & corrL, variable & corrK_cK, variable & corrL_cL) const
{
	if( !BuildFlux(fKL) ) return;

	get_variable<variable> Func(param);
	//advection part
	if( tag_U.isValid() )
	{
		//upwind data
		real       coefU  = fKL.Real            (tag_CONV_CU);
		real       rhsU   = fKL.Real            (tag_CONV_RU);
		Cell       cU     = Cell(m,fKL.Reference(tag_CONV_EU));
		//corrector in back cell
		real_array coefsK = fKL->RealArray      (tag_CONV_CB);
		ref_array  elemsK = fKL->ReferenceArray (tag_CONV_EB);
		real       rhsK   = fKL->Real           (tag_CONV_RB);
		//corrector in front cell
		real_array coefsL = fKL->RealArray      (tag_CONV_CF);
		ref_array  elemsL = fKL->ReferenceArray (tag_CONV_EF);
		real       rhsL   = fKL->Real           (tag_CONV_RF);
		//flag
		bulk       flag   = fKL->Bulk           (tag_CONV_F);

		fluxT += rhsU;
		if( cU.isValid() ) fluxT += Func(cU)*coefU;

		if( flag == 2 ) //both fluxes present
		{
			if( perform_correction_convection )
			{
				corrK += rhsK;
				corrL += rhsL;
				for(enumerator k = 0; k < elemsK.size(); ++k)
					if( elemsK[k].isValid() ) corrK += Func(elemsK[k])*coefsK[k];
				for(enumerator k = 0; k < elemsL.size(); ++k)
					if( elemsL[k].isValid() ) corrL += Func(elemsL[k])*coefsL[k];

				//if( scheme_type == NTPFA || scheme_type == NTPFA_PICARD )
				{
					if( elemsK.back().isValid() ) corrK_cK += Func(elemsK.back())*coefsK.back();
					if( elemsL.back().isValid() ) corrL_cL += Func(elemsL.back())*coefsL.back();
				}
			}
		}
		else
		{
			if( flag == 1 && perform_correction_convection ) // only corrK is present
			{
				corrK += rhsK;
				for(enumerator k = 0; k < elemsK.size(); ++k)
					if( elemsK[k].isValid() ) corrK += Func(elemsK[k])*coefsK[k];
			} //no corrections otherwise
		}
		
	} //advection part
}

int main(int argc, char ** argv)
{
	double nx = 2.0/7.0, ny = 6.0/7.0, nz = 3.0/7.0;
	double px = 0.5, py = 0.5, pz = 0.5;

	if( argc < 2 )
	{
		std::cout << "Usage: " << argv[0] << " mesh [mesh_out=grid.pmf] [nx=0] [ny=0] [nz=1] [px=0.5] [py=0.5] [pz=0.5]" << std::endl;
		return -1;
	}

	std::string grid_out = "grid.pmf";

	if( argc > 2 ) grid_out = std::string(argv[2]);
	if( argc > 3 ) nx = atof(argv[3]);
	if( argc > 4 ) ny = atof(argv[4]);
	if( argc > 5 ) nz = atof(argv[5]);
	if( argc > 6 ) px = atof(argv[6]);
	if( argc > 7 ) py = atof(argv[7]);
	if( argc > 8 ) pz = atof(argv[8]);

	double d = nx*px+ny*py+nz*pz;

	Mesh m;
	m.Load(argv[1]);
	m.SetTopologyCheck(NEED_TEST_CLOSURE|PROHIBIT_MULTILINE|PROHIBIT_MULTIPOLYGON|GRID_CONFORMITY|DEGENERATE_EDGE|DEGENERATE_FACE|DEGENERATE_CELL | FACE_EDGES_ORDER);
	//m.RemTopologyCheck(THROW_EXCEPTION);
	Tag sliced = m.CreateTag("SLICED",DATA_BULK,FACE|EDGE|NODE,FACE|EDGE|NODE,1);

	std::cout << "Cells: " << m.NumberOfCells() << std::endl;
	std::cout << "Faces: " << m.NumberOfFaces() << std::endl;

	MarkerType slice = m.CreateMarker();
	int nslice = 0, nmark = 0;
	for(Mesh::iteratorEdge it = m.BeginEdge(); it != m.EndEdge(); ++it)
	{
		double p[3];
		Storage::real_array c0 = it->getBeg()->Coords();
		Storage::real_array c1 = it->getEnd()->Coords();
		double r0 = c0[0]*nx+c0[1]*ny+c0[2]*nz - d;
		double r1 = c1[0]*nx+c1[1]*ny+c1[2]*nz - d;
		//std::cout << "r0 " << r0 << " r1 " << r1 << std::endl;
		if( r0*r1 < -1.0e-12 )
		{
			p[0] = (r0*c1[0] - r1*c0[0])/(r0-r1);
			p[1] = (r0*c1[1] - r1*c0[1])/(r0-r1);
			p[2] = (r0*c1[2] - r1*c0[2])/(r0-r1);
			//std::cout << "p " << p[0] << " " << p[1] << " " << p[2] << std::endl;
			Node n = m.CreateNode(p);
			n.Bulk(sliced) = 1;
			n.SetMarker(slice);
			bool was_sliced = it->HaveData(sliced) ? true : false;
			ElementArray<Edge> ret = Edge::SplitEdge(it->self(),ElementArray<Node>(&m,1,n.GetHandle()),0);
			if( was_sliced ) for(int q = 0; q < ret.size(); ++q) ret[q]->Bulk(sliced) = 1;
			nslice++;
		}
		else
		{
			if( fabs(r0) < 1.0e-6 )
			{
				it->getBeg()->SetMarker(slice);
				nmark++;
			}
			if( fabs(r1) < 1.0e-6 )
			{
				it->getEnd()->SetMarker(slice);
				nmark++;
			}
		}
	}

	std::cout << "sliced edges: " << nslice << " marked nodes: " << nmark << std::endl;

	if( !Element::CheckConnectivity(&m) )
		std::cout << "Connectivity is broken" << std::endl;
	
	nslice = 0;
	for(Mesh::iteratorFace it = m.BeginFace(); it != m.EndFace(); ++it)
	{
		ElementArray<Node> nodes = it->getNodes(slice); //those nodes should be ordered so that each pair forms an edge
		if( nodes.size() > 1 ) // if there is 1, then only one vertex touches the plane
		{
			//if there is more then two, then original face is non-convex
			if( nodes.size() > 2 ) std::cout << "Looks like face " << it->LocalID() << " is nonconvex" << std::endl;
			else
			{
				Edge e = m.CreateEdge(nodes).first;
				e.Bulk(sliced) = 1;
				e.SetMarker(slice);
				bool was_sliced = it->HaveData(sliced) ? true : false;
				ElementArray<Face> ret = Face::SplitFace(it->self(),ElementArray<Edge>(&m,1,e.GetHandle()),0);
				if( was_sliced ) for(int q = 0; q < ret.size(); ++q) ret[q]->Bulk(sliced) = 1;
				nslice++;
			}
		}
		//else std::cout << "Only one adjacent slice node, face " << it->LocalID() << std::endl;
	}

	nmark = 0;
	for(Mesh::iteratorEdge it = m.BeginEdge(); it != m.EndEdge(); ++it)
//.........这里部分代码省略.........

//.........这里部分代码省略.........
					T = lambdaK / dK;
				}
				else //un-splitted diffusion
				{
					uK = KKn;
					T = 0;
				}
						
				real bcconds[3] = {0.0,1.0,0.0}; //pure neumann boundary condition
				if( tag_BC.isValid() && fKL.HaveData(tag_BC) ) //are there boundary conditions on face?
				{
					//retrive boundary conditions
					real_array bc = fKL.RealArray(tag_BC);
					bcconds[0] = bc[0];
					bcconds[1] = bc[1];
					bcconds[2] = bc[2];
				}

				//account for boundary conditions
				R  = A*T*bcconds[2]/(bcconds[0] + bcconds[1]*T + 1.0e-100);
				uK *=A*  bcconds[0]/(bcconds[0] + bcconds[1]*T + 1.0e-100);
				T  *=A*  bcconds[0]/(bcconds[0] + bcconds[1]*T + 1.0e-100);
						
				//on BC
				if( uK.FrobeniusNorm() > 0 )
					flag_DIFF = 1;
				else
					flag_DIFF = 0;
			}
			else std::cout << "No adjacent cell on non-boundary face" << std::endl;
			//record data for diffusion part
			if( tag_K.isValid() )
			{
				fKL.Bulk(tag_DIFF_F) = flag_DIFF;
				fKL.Real(tag_DIFF_RT) = R;
				fKL.Real(tag_DIFF_CT) = T;
				real_array VT = fKL.RealArray(tag_DIFF_VT);
				VT[0] = uK(0,0);
				VT[1] = uK(0,1);
				VT[2] = uK(0,2);
				VT[3] = -uL(0,0);
				VT[4] = -uL(0,1);
				VT[5] = -uL(0,2);
			}
			//Advection part
			uK.Zero();
			uL.Zero();
			C = 0.0;
			cU = InvalidCell();
			R = 0;
			
			if( U > 0.0 ) //flow out of back cell to front cell
			{
				cU = cK;
				C = U*A;
				
				//upstream corrector
				uK = (xK - xKL)*U*A;
						
				//downstream corrector
				if( cL.isValid() ) //internal face
				{
					r0 = r = xK - xL;
					if( tag_K.isValid() ) //heterogeneous media
					{
						KD = KL - KK;