C++ Triangulation::numbOfTriangles方法代码示例

vector<double> FEM::computeDiscreteSolution(DiffusionReactionEqn const & PDE, 
                                            Triangulation& Omega,
                                            BoundaryConditions& BCs) {
	// data structures for final linear system A.xi = b:
	SymmetricCSlRMatrix A(Omega.generateAdjList()); // build final matrix portrait
	vector<double> b(Omega.numbOfNodes(), 0.), // load vector
	               xi(Omega.numbOfNodes(), 0.); // discrete solution	
	// data structures for assemby of A and b:
	SymmetricContainer<double> localMassMatrix(3), // for hat functions on triangles 
	                           localStiffnessMatrix(3), // we have 3 × 3 element matricies
	                           localRobinMatrix(2); // and 2 × 2 element matricies for Robin BCs (just like element matrix in 1D)
	array<double, 3> localLoadVector; // and their
	array<double, 2> localRobinVector; // friends, element vectors
	array<Node, 3> elementNodes, // nodes of the current triangle
	               elementMiddleNodes; // and nodes on the middle of edges
	array<Node, 2> edgeNodes; // nodes spanning an edge of the current triangle that is part of bndry
	Node midPoint; // middle point of the edge (to define which BCs to apply)
	double measure; // area of ith triangle / length of bndry edge of ith thiangle
	array<size_t, 3> l2g_elem; // local to global mapping of nodes on the element
	array<size_t, 2> l2g_edge; // and on the edge
	LocalIndex j, k, leftNodeIndex, rightNodeIndex; // dummy indicies
	for (size_t i = 0; i < Omega.numbOfTriangles(); ++i) {
		// (1) quadratures over elements
		// in order to assemble stiffness matrix and load vector,
		// it is convenient to iterate over mesh elements (i.e. triangles)
		elementNodes = Omega.getNodes(i); // get nodes of ith triangle
		for (j = 0; j < 3; ++j) // and middle nodes of its edges
			elementMiddleNodes[j] = elementNodes[k = nextIndex(j)].midPoint(elementNodes[nextIndex(k)]);
		measure = Omega.area(i); // compute area of ith triangle
		l2g_elem = Omega.l2g(i); // local to global mapping of nodes of ith element
		// compute
		// (1.1) local mass matrix,
		// (1.2) local stiffness matrix, and
		// (1.3) local load vector
		localStiffnessMatrix = computeLocalStiffnessMatrix(PDE.diffusionTerm(), elementNodes, elementMiddleNodes, measure);
		localMassMatrix      = computeLocalMassMatrix(PDE.reactionTerm(), elementNodes, measure);
		localLoadVector      = computeLocalLoadVector(PDE.forceTerm(), elementNodes, elementMiddleNodes, measure);
		// (1.4) assemble contributions
		for (j = 0; j < 3; ++j) {
			for (k = j; k < 3; ++k)
				A(l2g_elem[j], l2g_elem[k]) += localMassMatrix(j, k) + localStiffnessMatrix(j, k);
			b[l2g_elem[j]] += localLoadVector[j];
		}
		// (2) quadratures over edges
		// iterate over list of local indicies of boundary nodes
		for (LocalIndex edgeIndex : Omega.getBoundaryIndicies(i)) {
			// if edgeIndex = 2, then the edge against second node of ith triangle
			// is part of the boundary
			// so we need to assemble BCs here
			leftNodeIndex = nextIndex(edgeIndex); // local indicies of nodes that
			rightNodeIndex = nextIndex(leftNodeIndex); // define the edge
			edgeNodes = { elementNodes[leftNodeIndex], elementNodes[rightNodeIndex] }; // and the nodes themselves
			l2g_edge[0] = l2g_elem[leftNodeIndex]; // local to global nodes
			l2g_edge[1] = l2g_elem[rightNodeIndex]; // numeration mapping 
			measure = Omega.length(i, edgeIndex);
			// define BCs to apply
			midPoint = edgeNodes[0].midPoint(edgeNodes[1]);
			BCs.defineBCsAt(midPoint);
			// compute
			// (2.1) local Robin matrix
			// (2.2) local Robin vector
			localRobinMatrix = computeLocalRobinMatrix(BCs, edgeNodes, measure);
			localRobinVector = computeLocalRobinVector(BCs, edgeNodes, measure);
			// (2.3) assemble contributions
			for (j = 0; j < 2; ++j) {
				for (k = j; k < 2; ++k)
					A(l2g_edge[j], l2g_edge[k]) += localRobinMatrix(j, k);
				b[l2g_edge[j]] += localRobinVector[j];
			}
		}
	}
	// now we are ready to compute xi, A.xi = b
	xi = CG(A, b, xi, 10e-70);
	return xi;
}