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

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
}

//.........这里部分代码省略.........
			m->PrepareGeometricData(table);
		}
		BARRIER
		if( m->GetProcessorRank() == 0 ) std::cout << "Prepare geometric data: " << Timer()-ttt << std::endl;
		
		{
			Automatizator aut;
			Automatizator::MakeCurrent(&aut);
			INMOST_DATA_ENUM_TYPE iphi = aut.RegisterTag(phi,CELL);
			aut.EnumerateEntries();

			// Set the indeces intervals for the matrix and vectors
			R.SetInterval(aut.GetFirstIndex(),aut.GetLastIndex());
			R.InitLocks();
			Update.SetInterval(aut.GetFirstIndex(),aut.GetLastIndex());
			
			dynamic_variable Phi(aut,iphi);
			// Solve \nabla \cdot \nabla phi = f equation
			//for( Mesh::iteratorFace face = m->BeginFace(); face != m->EndFace(); ++face )
#if defined(USE_OMP)
#pragma omp parallel
#endif
			{
				variable flux; //should be more efficient to define here to avoid multiple memory allocations if storage for variations should be expanded
				rMatrix x1(3,1), x2(3,1), xf(3,1), n(3,1);
				double d1, d2, k1, k2, area, T, a, b, c;
#if defined(USE_OMP)
#pragma omp for
#endif
				for(Storage::integer iface = 0; iface < m->FaceLastLocalID(); ++iface ) if( m->isValidFace(iface) )
				{
					Face face = Face(m,ComposeFaceHandle(iface));
					Element::Status s1,s2;
					Cell r1 = face->BackCell();
					Cell r2 = face->FrontCell();
					if( ((!r1->isValid() || (s1 = r1->GetStatus()) == Element::Ghost)?0:1) +
						((!r2->isValid() || (s2 = r2->GetStatus()) == Element::Ghost)?0:1) == 0) continue;
					
					
					area = face->Area(); // Get the face area
					face->UnitNormal(n.data()); // Get the face normal
					face->Centroid(xf.data()); // Get the barycenter of the face
					r1->Centroid(x1.data());  // Get the barycenter of the cell
					k1 = n.DotProduct(rMatrix::FromTensor(tag_K[r1].data(),
														  tag_K[r1].size(),3)*n);
					d1 = fabs(n.DotProduct(xf-x1));
					if( !r2->isValid() ) // boundary condition
					{
						// bnd_pnt is a projection of the cell center to the face
						// a*pb + bT(pb-p1) = c
						// F = T(pb-p1)
						// pb = (c + bTp1)/(a+bT)
						// F = T/(a+bT)(c - ap1)
						T = k1/d1;
						a = 0;
						b = 1;
						c = 0;
						if( tag_BC.isValid() && face.HaveData(tag_BC) )
						{
							a = tag_BC[face][0];
							b = tag_BC[face][1];
							c = tag_BC[face][2];
							//std::cout << "a " << a << " b " << b << " c " << c << std::endl;
						}
						R.Lock(Phi.Index(r1));
						R[Phi.Index(r1)] -=  T/(a + b*T) * area * (c - a*Phi(r1));

ConvectionDiffusion::ConvectionDiffusion(Mesh * _m, Tag _tag_U, Tag _tag_K, Tag _tag_BC, MarkerType _boundary_face, bool correct_convection, bool correct_diffusion) 
		: m(_m), tag_U(_tag_U), tag_K(_tag_K), tag_BC(_tag_BC), boundary_face(_boundary_face)
{
	perform_correction_diffusion = correct_convection;
	perform_correction_convection = correct_diffusion;
	if( tag_U.isValid() )
	{
		//4 entries - triplet for upwind corrector, including the cell
		//back cell corrector
		tag_CONV_CB = m->CreateTag("CONV_BACK_COEFS"   ,DATA_REAL     ,FACE,NONE,4);
		tag_CONV_EB = m->CreateTag("CONV_BACK_ELEMS"   ,DATA_REFERENCE,FACE,NONE,4);
		tag_CONV_RB = m->CreateTag("CONV_BACK_RHS"     ,DATA_REAL     ,FACE,NONE,1);
		//front cell corrector
		tag_CONV_CF = m->CreateTag("CONV_FRONT_COEFS"  ,DATA_REAL     ,FACE,NONE,4);
		tag_CONV_EF = m->CreateTag("CONV_FRONT_ELEMS"  ,DATA_REFERENCE,FACE,NONE,4);
		tag_CONV_RF = m->CreateTag("CONV_FRONT_RHS"    ,DATA_REAL     ,FACE,NONE,1);
		//upstream
		tag_CONV_CU = m->CreateTag("CONV_UPSTREAM_COEF",DATA_REAL     ,FACE,NONE,1);
		tag_CONV_EU = m->CreateTag("CONV_UPSTREAM_ELEM",DATA_REFERENCE,FACE,NONE,1);
		tag_CONV_RU = m->CreateTag("CONV_UPSTREAM_RHS" ,DATA_REAL     ,FACE,NONE,1);
		tag_CONV_VU = m->CreateTag("CONV_UPSTREAM_CORR",DATA_REAL     ,FACE,NONE,6);
		//flag
		tag_CONV_F = m->CreateTag("CONV_FLAG"          ,DATA_BULK     ,FACE,NONE,1);
	}
	if( tag_K.isValid() )
	{
		//4 entries - triplet for transversal corrector, including the cell
		//back cell corrector
		tag_DIFF_CB = m->CreateTag("DIFF_BACK_COEFS"   ,DATA_REAL     ,FACE,NONE,4);
		tag_DIFF_EB = m->CreateTag("DIFF_BACK_ELEMS"   ,DATA_REFERENCE,FACE,NONE,4);
		tag_DIFF_RB = m->CreateTag("DIFF_BACK_RHS"     ,DATA_REAL     ,FACE,NONE,1);
		//front cell corrector
		tag_DIFF_CF = m->CreateTag("DIFF_FRONT_COEFS"  ,DATA_REAL     ,FACE,NONE,4);
		tag_DIFF_EF = m->CreateTag("DIFF_FRONT_ELEMS"  ,DATA_REFERENCE,FACE,NONE,4);
		tag_DIFF_RF = m->CreateTag("DIFF_FRONT_RHS"    ,DATA_REAL     ,FACE,NONE,1);
		//two-point transmissibility
		tag_DIFF_CT = m->CreateTag("DIFF_TP_COEF"      ,DATA_REAL     ,FACE,NONE,1);
		tag_DIFF_RT = m->CreateTag("DIFF_TP_RHS"       ,DATA_REAL     ,FACE,NONE,1);
		tag_DIFF_VT = m->CreateTag("DIFF_TP_CORR"      ,DATA_REAL     ,FACE,NONE,6);
		//flag
		tag_DIFF_F = m->CreateTag("DIFF_FLAG"          ,DATA_BULK     ,FACE,NONE,1);
	}

	build_flux = 0;
		
	if( m->GetProcessorsNumber() > 1 )
	{
		build_flux = m->CreateMarker();
#if defined(USE_OMP)
#pragma omp for
#endif
		for(integer q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
		{
			Face fKL = m->FaceByLocalID(q);
			Element::Status stat = fKL->BackCell()->GetStatus();
			if( fKL->FrontCell().isValid() ) stat |= fKL->FrontCell()->GetStatus();
			if( stat & (Element::Shared | Element::Owned) )
				fKL->SetMarker(build_flux);
		}	
	}

	initialization_failure = false; //there is a failure during intialization
	//Precompute two-point part, transversal correction direction,
	// upstream and upstream correction directions and flags
#if defined(USE_OMP)
#pragma omp parallel
#endif
	{
		//private memory
		rMatrix	xK(1,3), //center of cell K
				yK(1,3), //projection of xK onto interface
				xL(1,3), //center of cell L
				yL(1,3), //projection of xL onto interface
				xKL(1,3), //center of face common to K and L
				nKL(1,3), //normal vector to face
				KKn(1,3), //co-normal at cell K
				KLn(1,3), //co-normal at cell L
				v(1,3), //vector for correction
				uK(1,3), //vector for upwind correction in back cell
				uL(1,3), //vector for upwind correction in front cell
				r(1,3), //path to reach upstream concentration from downstream concentration
				r0(1,3),
				KK(3,3), //tenosr at cell K
				KL(3,3), //tensor at cell L
				KD(3,3), //difference of tensors
				gammaK(1,3), //transversal part of co-normal at cell K
				gammaL(1,3),  //transversal part of co-normal at cell L
				iT(3,3), //heterogeneous interpolation tensor
				iC(1,3) //heterogeneous interpolation correction
				;
		real	A, //area of the face
				U, //normal component of the velocity
				C, //coefficient for upstream cell
				T, //two-point transmissibility
				R, //right hand side
				dK, //distance from center to interface at cell K
				dL, //distance from center to interface at cell L
				lambdaK, //projection of co-normal onto normal at cell K
				lambdaL //projection of co-normal onto normal at cell L
				;
//.........这里部分代码省略.........