system_projection.c

一个用来实现偏微分方程中网格的计算库

                          // xz derivative
                          Ue(current_dof) = (gxplus(2) - gxminus(2))
                                            / 2. / TOLERANCE;
                          dof_is_fixed[current_dof] = true;
                          current_dof++;
                          // We need new points for yz
                          Point nyminus = elem->point(n),
                                nyplus = elem->point(n);
                          nyminus(1) -= TOLERANCE;
                          nyplus(1) += TOLERANCE;
                          Gradient gyminus = gptr(nyminus,
                                                  parameters,
                                                  system.name(),
                                                  system.variable_name(var));
                          Gradient gyplus = gptr(nyminus,
                                                 parameters,
                                                 system.name(),
                                                 system.variable_name(var));
                          // xz derivative
                          Ue(current_dof) = (gyplus(2) - gyminus(2))
                                            / 2. / TOLERANCE;
                          dof_is_fixed[current_dof] = true;
                          current_dof++;
                          // Getting a 2nd order xyz is more tedious
                          Point nxmym = elem->point(n),
                                nxmyp = elem->point(n),
                                nxpym = elem->point(n),
                                nxpyp = elem->point(n);
                          nxmym(0) -= TOLERANCE;
                          nxmym(1) -= TOLERANCE;
                          nxmyp(0) -= TOLERANCE;
                          nxmyp(1) += TOLERANCE;
                          nxpym(0) += TOLERANCE;
                          nxpym(1) -= TOLERANCE;
                          nxpyp(0) += TOLERANCE;
                          nxpyp(1) += TOLERANCE;
                          Gradient gxmym = gptr(nxmym,
                                                parameters,
                                                system.name(),
                                                system.variable_name(var));
                          Gradient gxmyp = gptr(nxmyp,
                                                parameters,
                                                system.name(),
                                                system.variable_name(var));
                          Gradient gxpym = gptr(nxpym,
                                                parameters,
                                                system.name(),
                                                system.variable_name(var));
                          Gradient gxpyp = gptr(nxpyp,
                                                parameters,
                                                system.name(),
                                                system.variable_name(var));
                          Number gxzplus = (gxpyp(2) - gxmyp(2))
                                         / 2. / TOLERANCE;
                          Number gxzminus = (gxpym(2) - gxmym(2))
                                          / 2. / TOLERANCE;
                          // xyz derivative
                          Ue(current_dof) = (gxzplus - gxzminus)
                                            / 2. / TOLERANCE;
                          dof_is_fixed[current_dof] = true;
                          current_dof++;
                        }
                    }
                }
              // Assume that other C_ONE elements have a single nodal
              // value shape function and nodal gradient component
              // shape functions
              else if (cont == C_ONE)
                {
                  libmesh_assert(nc == 1 + dim);
		  Ue(current_dof) = fptr(elem->point(n),
                                         parameters,
                                         system.name(),
                                         system.variable_name(var));
                  dof_is_fixed[current_dof] = true;
                  current_dof++;
                  Gradient g = gptr(elem->point(n),
                                    parameters,
                                    system.name(),
                                    system.variable_name(var));
                  for (unsigned int i=0; i!= dim; ++i)
                    {
		      Ue(current_dof) = g(i);
                      dof_is_fixed[current_dof] = true;
                      current_dof++;
                    }
                }
              else
                libmesh_error();
            }

          // In 3D, project any edge values next
          if (dim > 2 && cont != DISCONTINUOUS)
            for (unsigned int e=0; e != elem->n_edges(); ++e)
              {
		FEInterface::dofs_on_edge(elem, dim, fe_type, e,
                                          side_dofs);

                // Some edge dofs are on nodes and already
                // fixed, others are free to calculate
                unsigned int free_dofs = 0;
                for (unsigned int i=0; i != side_dofs.size(); ++i)
                  if (!dof_is_fixed[side_dofs[i]])
                    free_dof[free_dofs++] = i;

                // There may be nothing to project
                if (!free_dofs)
                  continue;

	        Ke.resize (free_dofs, free_dofs); Ke.zero();
	        Fe.resize (free_dofs); Fe.zero();
                // The new edge coefficients
                DenseVector<Number> Uedge(free_dofs);

                // Initialize FE data on the edge
                fe->attach_quadrature_rule (qedgerule.get());
	        fe->edge_reinit (elem, e);
	        const unsigned int n_qp = qedgerule->n_points();

	        // Loop over the quadrature points
	        for (unsigned int qp=0; qp<n_qp; qp++)
	          {
	            // solution at the quadrature point	      
	            Number fineval = fptr(xyz_values[qp],
                                          parameters,
                                          system.name(),
                                          system.variable_name(var));
	            // solution grad at the quadrature point	      
	            Gradient finegrad;
                    if (cont == C_ONE)
                      finegrad = gptr(xyz_values[qp], parameters,
                                      system.name(),
				      system.variable_name(var));

                    // Form edge projection matrix
                    for (unsigned int sidei=0, freei=0; 
                         sidei != side_dofs.size(); ++sidei)
                      {
                        unsigned int i = side_dofs[sidei];
                        // fixed DoFs aren't test functions
                        if (dof_is_fixed[i])
                          continue;
			for (unsigned int sidej=0, freej=0;
                             sidej != side_dofs.size(); ++sidej)
                          {
                            unsigned int j = side_dofs[sidej];
                            if (dof_is_fixed[j])
			      Fe(freei) -= phi[i][qp] * phi[j][qp] *
                                           JxW[qp] * Ue(j);
                            else
                              Ke(freei,freej) += phi[i][qp] *
						 phi[j][qp] * JxW[qp];
                            if (cont == C_ONE)
                              {
                                if (dof_is_fixed[j])
                                  Fe(freei) -= ((*dphi)[i][qp] *
					        (*dphi)[j][qp]) *
                                                JxW[qp] * Ue(j);
                                else
                                  Ke(freei,freej) += ((*dphi)[i][qp] *
                                                      (*dphi)[j][qp])
                                                      * JxW[qp];
                              }
                            if (!dof_is_fixed[j])
                              freej++;
                          }
                        Fe(freei) += phi[i][qp] * fineval * JxW[qp];
                        if (cont == C_ONE)
			  Fe(freei) += (finegrad * (*dphi)[i][qp]) *
                                       JxW[qp];
                        freei++;
                      }
	          }

                Ke.cholesky_solve(Fe, Uedge);

                // Transfer new edge solutions to element
		for (unsigned int i=0; i != free_dofs; ++i)
                  {
                    Number &ui = Ue(side_dofs[free_dof[i]]);
                    libmesh_assert(std::abs(ui) < TOLERANCE ||
                           std::abs(ui - Uedge(i)) < TOLERANCE);
                    ui = Uedge(i);
                    dof_is_fixed[side_dofs[free_dof[i]]] = true;
                  }
              }
                 
	  // Project any side values (edges in 2D, faces in 3D)
          if (dim > 1 && cont != DISCONTINUOUS)
            for (unsigned int s=0; s != elem->n_sides(); ++s)
              {
		FEInterface::dofs_on_side(elem, dim, fe_type, s,
                                          side_dofs);

		// Some side dofs are on nodes/edges and already
                // fixed, others are free to calculate
                unsigned int free_dofs = 0;
                for (unsigned int i=0; i != side_dofs.size(); ++i)
                  if (!dof_is_fixed[side_dofs[i]])
                    free_dof[free_dofs++] = i;

                // There may be nothing to project
                if (!free_dofs)
                  continue;
             
	        Ke.resize (free_dofs, free_dofs); Ke.zero();
	        Fe.resize (free_dofs); Fe.zero();
                // The new side coefficients
                DenseVector<Number> Uside(free_dofs);

                // Initialize FE data on the side
                fe->attach_quadrature_rule (qsiderule.get());
	        fe->reinit (elem, s);
	        const unsigned int n_qp = qsiderule->n_points();

	        // Loop over the quadrature points
	        for (unsigned int qp=0; qp<n_qp; qp++)
	          {
	            // solution at the quadrature point	      
	            Number fineval = fptr(xyz_values[qp],
                                          parameters,
                                          system.name(),
                                          system.variable_name(var));
	            // solution grad at the quadrature point	      
	            Gradient finegrad;
                    if (cont == C_ONE)
                      finegrad = gptr(xyz_values[qp], parameters,
                                      system.name(),
				      system.variable_name(var));

                    // Form side projection matrix
                    for (unsigned int sidei=0, freei=0;
                         sidei != side_dofs.size(); ++sidei)
                      {
                        unsigned int i = side_dofs[sidei];
                        // fixed DoFs aren't test functions
                        if (dof_is_fixed[i])
                          continue;
			for (unsigned int sidej=0, freej=0;
                             sidej != side_dofs.size(); ++sidej)
                          {
                            unsigned int j = side_dofs[sidej];
                            if (dof_is_fixed[j])
			      Fe(freei) -= phi[i][qp] * phi[j][qp] *
                                           JxW[qp] * Ue(j);
                            else
			      Ke(freei,freej) += phi[i][qp] *
						 phi[j][qp] * JxW[qp];
                            if (cont == C_ONE)
                              {
                                if (dof_is_fixed[j])
                                  Fe(freei) -= ((*dphi)[i][qp] *
					        (*dphi)[j][qp]) *
                                               JxW[qp] * Ue(j);
                                else
                                  Ke(freei,freej) += ((*dphi)[i][qp] *
						      (*dphi)[j][qp])
                                                     * JxW[qp];
                              }
                            if (!dof_is_fixed[j])
                              freej++;
                          }
                        Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
                        if (cont == C_ONE)
			  Fe(freei) += (finegrad * (*dphi)[i][qp]) *
                                       JxW[qp];
                        freei++;
                      }
	          }

                Ke.cholesky_solve(Fe, Uside);

                // Transfer new side solutions to element
		for (unsigned int i=0; i != free_dofs; ++i)
                  {
                    Number &ui = Ue(side_dofs[free_dof[i]]);
                    libmesh_assert(std::abs(ui) < TOLERANCE ||
                           std::abs(ui - Uside(i)) < TOLERANCE);
                    ui = Uside(i);
                    dof_is_fixed[side_dofs[free_dof[i]]] = true;
                  }
              }

	  // Project the interior values, finally

	  // Some interior dofs are on nodes/edges/sides and
          // already fixed, others are free to calculate
          unsigned int free_dofs = 0;
          for (unsigned int i=0; i != n_dofs; ++i)
            if (!dof_is_fixed[i])
              free_dof[free_dofs++] = i;

          // There may be nothing to project
          if (free_dofs)
            {
             
	  Ke.resize (free_dofs, free_dofs); Ke.zero();
	  Fe.resize (free_dofs); Fe.zero();
          // The new interior coefficients
          DenseVector<Number> Uint(free_dofs);

          // Initialize FE data
          fe->attach_quadrature_rule (qrule.get());
	  fe->reinit (elem);
	  const unsigned int n_qp = qrule->n_points();

	  // Loop over the quadrature points
	  for (unsigned int qp=0; qp<n_qp; qp++)
	    {
	      // solution at the quadrature point	      
	      Number fineval = fptr(xyz_values[qp],
                                    parameters,
                                    system.name(),
                                    system.variable_name(var));
	      // solution grad at the quadrature point	      
	      Gradient finegrad;
              if (cont == C_ONE)
                finegrad = gptr(xyz_values[qp], parameters,
                                system.name(),
				system.variable_name(var));

              // Form interior projection matrix
              for (unsigned int i=0, freei=0; i != n_dofs; ++i)
                {
                  // fixed DoFs aren't test functions
                  if (dof_is_fixed[i])
                    continue;
		  for (unsigned int j=0, freej=0; j != n_dofs; ++j)
                    {
                      if (dof_is_fixed[j])
			Fe(freei) -= phi[i][qp] * phi[j][qp] * JxW[qp]
                                     * Ue(j);
                      else
			Ke(freei,freej) += phi[i][qp] * phi[j][qp] *
                                           JxW[qp];
                      if (cont == C_ONE)
                        {
                          if (dof_is_fixed[j])
			    Fe(freei) -= ((*dphi)[i][qp] *
					 (*dphi)[j][qp]) * JxW[qp] *
                                         Ue(j);
                          else
			    Ke(freei,freej) += ((*dphi)[i][qp] *
						(*dphi)[j][qp]) *
                                               JxW[qp];
                        }
                      if (!dof_is_fixed[j])
                        freej++;
                    }
		  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
                  if (cont == C_ONE)
                    Fe(freei) += (finegrad * (*dphi)[i][qp]) * JxW[qp];
                  freei++;
                }
	    }
          Ke.cholesky_solve(Fe, Uint);

          // Transfer new interior solutions to element
	  for (unsigned int i=0; i != free_dofs; ++i)
            {
              Number &ui = Ue(free_dof[i]);
              libmesh_assert(std::abs(ui) < TOLERANCE ||
                     std::abs(ui - Uint(i)) < TOLERANCE);
              ui = Uint(i);
              dof_is_fixed[free_dof[i]] = true;
            }

            } // if there are free interior dofs

          // Make sure every DoF got reached!
	  for (unsigned int i=0; i != n_dofs; ++i)
            libmesh_assert(dof_is_fixed[i]);

	  const unsigned int
	    first = new_vector.first_local_index(),
	    last  = new_vector.last_local_index();
	  
	  // Lock the new_vector since it is shared among threads.
	  {
	    Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
	    
	    for (unsigned int i = 0; i < n_dofs; i++) 
	      // We may be projecting a new zero value onto
	      // an old nonzero approximation - RHS
	      // if (Ue(i) != 0.)
	      if ((dof_indices[i] >= first) &&
		  (dof_indices[i] <  last))
                {
                  new_vector.set(dof_indices[i], Ue(i));
                }
	  }
        }  // end elem loop
    } // end variables loop
}