Commented assembly as implemented by FreeFem++

MatriceMorse_Build.cpp

template <class R>
template <class FESpace> 
void MatriceMorse<R>::Build(const FESpace & Uh, const FESpace & Vh, bool sym)
{
	typedef typename FESpace::Mesh Mesh;

	symetrique = sym;

	// Coefficients
	a = 0;

	// Lines
	lg = 0;

	// Columns
	cl = 0;

  	// same Mesh
	ffassert( &Uh.Th == &Vh.Th);
	const Mesh & Th(Uh.Th);

	int nbe = Uh.NbOfElements;
	int nbn_u = Uh.NbOfNodes;
	int nbn_v = Vh.NbOfNodes;

	KN<int> pe_u(nbn_u + 1 + Uh.SizeToStoreAllNodeofElement());

	// Reverse connectivity table encoded as a single array:
	// * first nbn_u values are start/end offsets in the second half of the pe_u
	// * second half of pe_u contains elements id for each vertex k: pe_u[pe_u[k]] to pe_u[pe_u[k+1]]

	// First, count #elements per node. A node can be a vertex, an edge, a triangle, a tetrahedron
	pe_u = 0;
	for (int k = 0; k < nbe; k++)
	{ 
		// Number of node for element k, for example if Uh(k) is a triangle 
		// 3 if it has 1 dof per verts
		// 6 if it has 1 dof per verts and edges ...
		int nbne = Uh(k);
		for (int in = 0; in < nbne; in++)
		{
			// Node Uh(k,in) belongs to pe_u[Uh(k,in) + 1] distinct elements
			pe_u[Uh(k,in) + 1]++;
		}
	}

	// Then "prefix sum" the count array: this gives offsets in the second half of pe_u
	// The offsets are actually shifted one cell to the left, it will be unshifted when the second half of pe_u is populated.
	int kk = nbn_u + 1, kkk = kk;
	pe_u[0] = kk;
	for (int in1 = 1; in1 <= nbn_u; in1++)
	{
		kk += pe_u[in1];
		pe_u[in1] = kkk;
		kkk = kk;
	} 


	ffassert(kk == nbn_u + 1 + Uh.SizeToStoreAllNodeofElement());

	// Finaly, we populate the second half of pe_u with elements id.
	// Incrementing pe_u[Uh(k,in) + 1] keeps track of the next cell used to store the next element id for vertex Uh(k,in)
	// as well as shifting back to the right the first half of pe_u so that we access Uh(k,in)'s elememts doing pe_u[pe_u[Uh(k,in)]
	// instead of pe_u[pe_u[Uh(k,in) + 1]
	for (int k = 0; k < nbe; k++)
	{ 
		int nbne = Uh(k);
		for (int in = 0; in < nbne; in++)
		{
			pe_u[pe_u[Uh(k,in) + 1]++] = k;
		}
	}


	// Now it will scan the lines of the linear system 2 times:
	// * the first pass (step == 0) it counts the number of non zero coefficients in order to allocate memory for arrays a and cl
	// * the second pass (step == 1) it populates the cl array with DoF ids

	// Mark and color are used to remenber that we already accounted for a non zero coeff (column) for the current line.
	// color is incremented every line
	KN<int> mark(nbn_v);
	int color = 0;
	mark = color++;

	lg = new int [this->n + 1];
	ffassert(lg);

	for (int step = 0; step < 2; step++) 
	{ 

		int ilg = 0;
		lg[0] = ilg;
		int kij = 0;

		// foreach line "in"
		for (int in = 0; in < nbn_u; in++)
		{

			// Number of non-zero coeff for that line
			int nbj = 0;

			// Offset in the cl array for that line
			int kijs = kij;

			// foreach element containing node "in" (using the reverse connectivity table)
			for (int kk = pe_u[in]; kk < pe_u[in+1]; kk++)
			{
				// ke is the id of an lement containing node "in"
				int ke = pe_u[kk];

				// BEGINING OF MISTERY CODE RELATED TO FINITE VOLUME METHOD
				int tabk[10];
				int ltab = 0;
				tabk[ltab++] = ke;
				if( VF )
				{
					ltab += Th.GetAllElementAdj(ke,tabk+ltab);
				}
				tabk[ltab]=-1;
				// END OF MISTERY CODE RELATED TO FINITE VOLUME METHOD

				
				// MISTERY LOOP: for the FINITE *ELEMENT* METHOD, 
				// ltab is equal to *one* hence the loop is executed *once*
				for(int ik = 0; ik < ltab; ++ik)
				{

					// for the FINITE *ELEMENT* METHOD k == tabk[0] is always equal to ke
					// so that we only count columns/dof related to element ke
					int k = tabk[ik];
					throwassert(k >= 0 && k < nbe);

					// foreach node of element k
					int njloc = Vh(k);
					for (int jloc = 0; jloc < njloc; jloc++)
					{ 
						// jn is the id of the other node
						int jn = Vh(k,jloc);

						// if jn hasn't been accounted for already
						// we should account for it only if the matrix is non symmetrical matrix or it's a coeff of the lower triangular block of the matrix
						// Question " why "jn < in" rather than "jn <= in" ?
						// Answer: not all DoF coupling for (in,in) are part of the lower triangular block of the matrix, we need to be extra careful
						if (mark[jn] != color && (!sym ||  jn < in)) 
						{
							// Rememeber we accounted for nj
							mark[jn] = color;

							// We account for degrees of freedom, not nodes
							// so that we must retrieve the #DoF for node jn
							int fdf=Vh.FirstDFOfNode(jn);
							int ldf=Vh.LastDFOfNode(jn);

							// During the second pass, cl has been allocated and can be populated
							if (step)
							{
								// j is a DoF id
								for (int j = fdf; j < ldf; j++)
								{
									// Increment kij to keep track of the next cell for the next non zero
									cl[kij++] = j;
								}
							}

							// Update # of non zero coeff for node "in" by adding the # of DoF for jn
							nbj += ldf-fdf;
						}
					} // END OF COLUMN LOOP
		    		} // END OF MISTERY LOOP
		  	} // END OF ELEMENT LOOP


			int fdf = Uh.FirstDFOfNode(in);
			int ldf = Uh.LastDFOfNode(in);
			int kijl = kij;

			if (step)
			{
				// During the second pass, at the end of the line loop
				// the non zero indices are sorted in increasing order
				HeapSort(cl+kijs,kij-kijs);

				// Account for non zero as many times as there are DoF related to vertex "in"
				for (int i = fdf; i < ldf; i++)
				{ 
					//  Copy the line if not the first as the first one has been populated in the COLUMN LOOP
					if (i != fdf)
					{	
						// Duplicate
						for (int k = kijs; k < kijl; k++)
						{
							cl[kij++] = cl[k];
						}
					}

					// We did not account for the diagonal couplings in a first place
					// See Question/Answer in the COLUMN loop
					if (sym)
					{
						for(int j = fdf; j <= i; j++)
						{
							cl[kij++] = j;
						}
					}

					throwassert(kij==lg[i+1]);// verif           
				}
			}
		  	else
		  	{
	  		    // Account for non zero as many times as there are DoF related to vertex "in"
			    for (int i = fdf; i < ldf; i++)
			    { 			
			    	// We did not account for the diagonal couplings in a first place
			    	// See Question/Answer in the COLUMN loop
		      		if (sym) 
					{						
			      		ilg += ++nbj;
			      	}
			      	else 
					{
		      			ilg += nbj;
			      	}
			      	
			      	lg[i+1] = ilg;
			    }
			}

			// Next color
			color++;

		} // END OF ELEMENT LOOP

		// At the end of pass 1, we know the number of non zero elements so that we can allocate memory for a and cl
		if (step == 0) 
		{
			// Alloc
			nbcoef = ilg;		 
			a = new R[nbcoef];
			cl = new int [nbcoef];
		
			ffassert(a && cl);

			// "Zero" a
			for (int i = 0; i < nbcoef; i++)
			{ 
			  a[i] = 0;
			}
		}
	} // END OF STEP LOOP
}