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

int main(int argc,char ** argv)
{
	Solver::Initialize(&argc,&argv,""); // Initialize the solver and MPI activity
#if defined(USE_PARTITIONER)
	Partitioner::Initialize(&argc,&argv); // Initialize the partitioner activity
#endif
	if( argc > 1 )
	{
		TagReal phi;
		TagReal tag_F;
		TagRealArray tag_K;
		TagRealArray tag_BC;
		TagReal phi_ref;
		Mesh * m = new Mesh(); // Create an empty mesh
		double ttt = Timer();
		bool repartition = false;
		m->SetCommunicator(INMOST_MPI_COMM_WORLD); // Set the MPI communicator for the mesh
		if( m->GetProcessorRank() == 0 ) // If the current process is the master one
			std::cout << argv[0] << std::endl;

		if( m->isParallelFileFormat(argv[1]) )
		{
			m->Load(argv[1]); // Load mesh from the parallel file format
			repartition = true;
		}
		else
		{
			if( m->GetProcessorRank() == 0 )
				m->Load(argv[1]); // Load mesh from the serial file format
		}
		BARRIER;
		if( m->GetProcessorRank() == 0 ) std::cout << "Processors: " << m->GetProcessorsNumber() << std::endl;
		if( m->GetProcessorRank() == 0 ) std::cout << "Load(MPI_File): " << Timer()-ttt << std::endl;

		//~ double ttt2 = Timer();
		//~ Mesh t;
		//~ t.SetCommunicator(INMOST_MPI_COMM_WORLD);
		//~ t.SetParallelFileStrategy(0);
		//~ t.Load(argv[1]);
		//~ BARRIER
		//~ if( m->GetProcessorRank() == 0 ) std::cout << "Load(MPI_Scatter): " << Timer()-ttt2 << std::endl;

#if defined(USE_PARTITIONER)
		if (m->GetProcessorsNumber() > 1)
		{ // currently only non-distributed meshes are supported by Inner_RCM partitioner
			ttt = Timer();
			Partitioner * p = new Partitioner(m);
			p->SetMethod(Partitioner::INNER_KMEANS,Partitioner::Partition); // Specify the partitioner
			p->Evaluate(); // Compute the partitioner and store new processor ID in the mesh
			delete p;
			BARRIER;

			if( m->GetProcessorRank() == 0 ) std::cout << "Evaluate: " << Timer()-ttt << std::endl;

			ttt = Timer();
			m->Redistribute(); // Redistribute the mesh data
			m->ReorderEmpty(CELL|FACE|EDGE|NODE); // Clean the data after reordring
			BARRIER;

			if( m->GetProcessorRank() == 0 ) std::cout << "Redistribute: " << Timer()-ttt << std::endl;
		}
#endif

		ttt = Timer();
		phi = m->CreateTag("Solution",DATA_REAL,CELL,NONE,1); // Create a new tag for the solution phi
		
		bool makerefsol = true;
		
		if( m->HaveTag("PERM" ) )
		{
			tag_K = m->GetTag("PERM");
			makerefsol = false;
			std::cout << "Permeability from grid" << std::endl;
		}
		else
		{
			std::cout << "Set perm" << std::endl;
			tag_K = m->CreateTag("PERM",DATA_REAL,CELL,NONE,1); // Create a new tag for K tensor
			for( Mesh::iteratorCell cell = m->BeginCell(); cell != m->EndCell(); ++cell ) // Loop over mesh cells
				tag_K[*cell][0] = 1.0; // Store the tensor K value into the tag
		}
		
		
		
		if( m->HaveTag("BOUNDARY_CONDITION") )
		{
			tag_BC = m->GetTag("BOUNDARY_CONDITION");
			makerefsol = false;
			std::cout << "Boundary conditions from grid" << std::endl;
		}
		else
		{
			std::cout << "Set boundary conditions" << std::endl;
			double x[3];
			tag_BC = m->CreateTag("BOUNDARY_CONDITION",DATA_REAL,FACE,FACE,3);
			for( Mesh::iteratorFace face = m->BeginFace(); face != m->EndFace(); ++face )
				if( face->Boundary() && !(face->GetStatus() == Element::Ghost) )
				{
					face->Centroid(x);
					tag_BC[*face][0] = 1; //dirichlet
//.........这里部分代码省略.........

//.........这里部分代码省略.........
				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
				;
		const real eps = degenerate_diffusion_regularization;
		Cell cK, cL, cU;
		Face fKL;
		bulk flag_DIFF, flag_CONV;
		KK.Zero();
		KL.Zero();
		KD.Zero();
		U = 0.0;
#if defined(USE_OMP)
#pragma omp for
#endif
		for(integer q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
		{
			fKL = m->FaceByLocalID(q);
			if( !BuildFlux(fKL) ) continue;

			fKL.Centroid(xKL.data());
			fKL.UnitNormal(nKL.data());
			A = fKL.Area();
			if( tag_U.isValid() ) U = fKL.Real(tag_U);
					
			cK = fKL.BackCell();
			cL = fKL.FrontCell();
			assert(cK.isValid());
					
			cK.Centroid(xK.data());
			dK = nKL.DotProduct(xKL-xK);
			yK = xK + dK*nKL;
			if( tag_K.isValid() )
				KK = rMatrix::FromTensor(cK.RealArray(tag_K).data(),cK.RealArray(tag_K).size());//.Transpose();
			KKn = nKL*KK;
			lambdaK = nKL.DotProduct(KKn);
			gammaK = KKn - lambdaK*nKL;
					
			//Diffusion part
			uK.Zero();
			uL.Zero();
			if( cL.isValid() ) //internal, both cells are present
			{
				cL.Centroid(xL.data());
				dL = nKL.DotProduct(xL-xKL);
				yL = xL - dL*nKL;
				if( tag_K.isValid() )
					KL = rMatrix::FromTensor(cL.RealArray(tag_K).data(),cL.RealArray(tag_K).size());//.Transpose();
				KLn = nKL*KL;
				lambdaL = nKL.DotProduct(KLn);
				gammaL = KLn - lambdaL*nKL;