libMesh::ProjectFEMSolution Class Reference

Public Member Functions

 ProjectFEMSolution (const System &system_in, FEMFunctionBase< Number > *f_in, FEMFunctionBase< Gradient > *g_in, NumericVector< Number > &new_v_in)
 
 ProjectFEMSolution (const ProjectFEMSolution &in)
 
void operator() (const ConstElemRange &range) const
 

Private Attributes

const Systemsystem
 
AutoPtr< FEMFunctionBase
< Number > > 
f
 
AutoPtr< FEMFunctionBase
< Gradient > > 
g
 
NumericVector< Number > & new_vector
 

Detailed Description

This class implements projecting an arbitrary function to the current mesh. This may be exectued in parallel on multiple threads.

Definition at line 160 of file system_projection.C.

Constructor & Destructor Documentation

libMesh::ProjectFEMSolution::ProjectFEMSolution ( const System system_in,
FEMFunctionBase< Number > *  f_in,
FEMFunctionBase< Gradient > *  g_in,
NumericVector< Number > &  new_v_in 
)
inline

Definition at line 170 of file system_projection.C.

References f, and libMesh::libmesh_assert().

173  :
174  system(system_in),
175  f(f_in ? f_in->clone() : AutoPtr<FEMFunctionBase<Number> >(NULL)),
176  g(g_in ? g_in->clone() : AutoPtr<FEMFunctionBase<Gradient> >(NULL)),
177  new_vector(new_v_in)
178  {
179  libmesh_assert(f.get());
180  }
libMesh::ProjectFEMSolution::ProjectFEMSolution ( const ProjectFEMSolution in)
inline

Definition at line 182 of file system_projection.C.

References f, and libMesh::libmesh_assert().

182  :
183  system(in.system),
184  f(in.f.get() ? in.f->clone() : AutoPtr<FEMFunctionBase<Number> >(NULL)),
185  g(in.g.get() ? in.g->clone() : AutoPtr<FEMFunctionBase<Gradient> >(NULL)),
186  new_vector(in.new_vector)
187  {
188  libmesh_assert(f.get());
189  }

Member Function Documentation

void libMesh::ProjectFEMSolution::operator() ( const ConstElemRange range) const

This method projects an arbitrary solution to the current mesh. The input function f gives the arbitrary solution, while the new_vector (which should already be correctly sized) gives the solution (to be computed) on the current mesh.

Definition at line 1833 of file system_projection.C.

References std::abs(), libMesh::Variable::active_on_subdomain(), libMesh::StoredRange< iterator_type, object_type >::begin(), libMeshEnums::C_ONE, libMeshEnums::C_ZERO, libMesh::DenseMatrix< T >::cholesky_solve(), libMesh::FunctionBase< Output >::component(), libMesh::dim, libMeshEnums::DISCONTINUOUS, libMesh::FEInterface::dofs_on_edge(), libMesh::FEInterface::dofs_on_side(), libMesh::StoredRange< iterator_type, object_type >::end(), libMesh::FEType::family, libMesh::NumericVector< T >::first_local_index(), libMesh::FEAbstract::get_continuity(), libMesh::System::get_dof_map(), libMesh::FEGenericBase< T >::get_dphi(), libMesh::FEAbstract::get_fe_type(), libMesh::FEAbstract::get_JxW(), libMesh::FEGenericBase< T >::get_phi(), libMesh::FEAbstract::get_xyz(), libMeshEnums::HERMITE, libMesh::Elem::is_vertex(), libMesh::NumericVector< T >::last_local_index(), libMesh::libmesh_assert(), libMesh::FEInterface::n_dofs_at_node(), libMesh::Elem::n_edges(), n_nodes, libMesh::Elem::n_nodes(), libMesh::QBase::n_points(), libMesh::Elem::n_sides(), libMesh::Elem::point(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMeshEnums::SCALAR, libMesh::NumericVector< T >::set(), libMesh::Threads::spin_mtx, libMesh::Elem::subdomain_id(), system, libMesh::System::time, libMesh::TOLERANCE, libMesh::Variable::type(), libMesh::Elem::type(), libMesh::DofMap::variable(), libMesh::System::variable_scalar_number(), libMesh::DenseMatrix< T >::zero(), and libMesh::DenseVector< T >::zero().

1834 {
1835  // We need data to project
1836  libmesh_assert(f.get());
1837 
1845  FEMContext context( system );
1846 
1847  // The number of variables in this system
1848  const unsigned int n_variables = context.n_vars();
1849 
1850  // The dimensionality of the current mesh
1851  const unsigned int dim = context.get_dim();
1852 
1853  // The DofMap for this system
1854  const DofMap& dof_map = system.get_dof_map();
1855 
1856  // The element matrix and RHS for projections.
1857  // Note that Ke is always real-valued, whereas
1858  // Fe may be complex valued if complex number
1859  // support is enabled
1860  DenseMatrix<Real> Ke;
1861  DenseVector<Number> Fe;
1862  // The new element coefficients
1863  DenseVector<Number> Ue;
1864 
1865  // FIXME: Need to generalize this to vector-valued elements. [PB]
1866  FEBase* fe = NULL;
1867  FEBase* side_fe = NULL;
1868  FEBase* edge_fe = NULL;
1869 
1870  // First, loop over all variables and make sure the shape functions phi will
1871  // get built as well as the shape function gradients if the gradient functor
1872  // is supplied.
1873  for (unsigned int var=0; var<n_variables; var++)
1874  {
1875  context.get_element_fe( var, fe );
1876  if (fe->get_fe_type().family == SCALAR)
1877  continue;
1878  if( dim > 1 )
1879  context.get_side_fe( var, side_fe );
1880  if( dim > 2 )
1881  context.get_edge_fe( var, edge_fe );
1882 
1883  fe->get_phi();
1884  if( dim > 1 )
1885  side_fe->get_phi();
1886  if( dim > 2 )
1887  edge_fe->get_phi();
1888 
1889  const FEContinuity cont = fe->get_continuity();
1890  if(cont == C_ONE)
1891  {
1892  libmesh_assert(g.get());
1893  fe->get_dphi();
1894  side_fe->get_dphi();
1895  if( dim > 2 )
1896  edge_fe->get_dphi();
1897  }
1898  }
1899 
1900  // Now initialize any user requested shape functions
1901  f->init_context(context);
1902  if(g.get())
1903  g->init_context(context);
1904 
1905  std::vector<unsigned int> side_dofs;
1906 
1907  // Iterate over all the elements in the range
1908  for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it != range.end(); ++elem_it)
1909  {
1910  const Elem* elem = *elem_it;
1911 
1912  context.pre_fe_reinit(system, elem);
1913 
1914  // Loop over all the variables in the system
1915  for (unsigned int var=0; var<n_variables; var++)
1916  {
1917  const Variable& variable = dof_map.variable(var);
1918 
1919  const FEType& fe_type = variable.type();
1920 
1921  if (fe_type.family == SCALAR)
1922  continue;
1923 
1924  // Per-subdomain variables don't need to be projected on
1925  // elements where they're not active
1926  if (!variable.active_on_subdomain(elem->subdomain_id()))
1927  continue;
1928 
1929  const FEContinuity cont = fe->get_continuity();
1930 
1931  const unsigned int var_component =
1933 
1934  const std::vector<dof_id_type>& dof_indices =
1935  context.get_dof_indices(var);
1936 
1937  // The number of DOFs on the element
1938  const unsigned int n_dofs =
1939  libmesh_cast_int<unsigned int>(dof_indices.size());
1940 
1941  // Fixed vs. free DoFs on edge/face projections
1942  std::vector<char> dof_is_fixed(n_dofs, false); // bools
1943  std::vector<int> free_dof(n_dofs, 0);
1944 
1945  // The element type
1946  const ElemType elem_type = elem->type();
1947 
1948  // The number of nodes on the new element
1949  const unsigned int n_nodes = elem->n_nodes();
1950 
1951  // Zero the interpolated values
1952  Ue.resize (n_dofs); Ue.zero();
1953 
1954  // In general, we need a series of
1955  // projections to ensure a unique and continuous
1956  // solution. We start by interpolating nodes, then
1957  // hold those fixed and project edges, then
1958  // hold those fixed and project faces, then
1959  // hold those fixed and project interiors
1960 
1961  // Interpolate node values first
1962  unsigned int current_dof = 0;
1963  for (unsigned int n=0; n!= n_nodes; ++n)
1964  {
1965  // FIXME: this should go through the DofMap,
1966  // not duplicate dof_indices code badly!
1967  const unsigned int nc =
1968  FEInterface::n_dofs_at_node (dim, fe_type, elem_type,
1969  n);
1970  if (!elem->is_vertex(n))
1971  {
1972  current_dof += nc;
1973  continue;
1974  }
1975  if (cont == DISCONTINUOUS)
1976  {
1977  libmesh_assert_equal_to (nc, 0);
1978  }
1979  // Assume that C_ZERO elements have a single nodal
1980  // value shape function
1981  else if (cont == C_ZERO)
1982  {
1983  libmesh_assert_equal_to (nc, 1);
1984  Ue(current_dof) = f->component(context,
1985  var_component,
1986  elem->point(n),
1987  system.time);
1988  dof_is_fixed[current_dof] = true;
1989  current_dof++;
1990  }
1991  // The hermite element vertex shape functions are weird
1992  else if (fe_type.family == HERMITE)
1993  {
1994  Ue(current_dof) = f->component(context,
1995  var_component,
1996  elem->point(n),
1997  system.time);
1998  dof_is_fixed[current_dof] = true;
1999  current_dof++;
2000  Gradient grad = g->component(context,
2001  var_component,
2002  elem->point(n),
2003  system.time);
2004  // x derivative
2005  Ue(current_dof) = grad(0);
2006  dof_is_fixed[current_dof] = true;
2007  current_dof++;
2008  if (dim > 1)
2009  {
2010  // We'll finite difference mixed derivatives
2011  Point nxminus = elem->point(n),
2012  nxplus = elem->point(n);
2013  nxminus(0) -= TOLERANCE;
2014  nxplus(0) += TOLERANCE;
2015  Gradient gxminus = g->component(context,
2016  var_component,
2017  nxminus,
2018  system.time);
2019  Gradient gxplus = g->component(context,
2020  var_component,
2021  nxplus,
2022  system.time);
2023  // y derivative
2024  Ue(current_dof) = grad(1);
2025  dof_is_fixed[current_dof] = true;
2026  current_dof++;
2027  // xy derivative
2028  Ue(current_dof) = (gxplus(1) - gxminus(1))
2029  / 2. / TOLERANCE;
2030  dof_is_fixed[current_dof] = true;
2031  current_dof++;
2032 
2033  if (dim > 2)
2034  {
2035  // z derivative
2036  Ue(current_dof) = grad(2);
2037  dof_is_fixed[current_dof] = true;
2038  current_dof++;
2039  // xz derivative
2040  Ue(current_dof) = (gxplus(2) - gxminus(2))
2041  / 2. / TOLERANCE;
2042  dof_is_fixed[current_dof] = true;
2043  current_dof++;
2044  // We need new points for yz
2045  Point nyminus = elem->point(n),
2046  nyplus = elem->point(n);
2047  nyminus(1) -= TOLERANCE;
2048  nyplus(1) += TOLERANCE;
2049  Gradient gyminus = g->component(context,
2050  var_component,
2051  nyminus,
2052  system.time);
2053  Gradient gyplus = g->component(context,
2054  var_component,
2055  nyplus,
2056  system.time);
2057  // xz derivative
2058  Ue(current_dof) = (gyplus(2) - gyminus(2))
2059  / 2. / TOLERANCE;
2060  dof_is_fixed[current_dof] = true;
2061  current_dof++;
2062  // Getting a 2nd order xyz is more tedious
2063  Point nxmym = elem->point(n),
2064  nxmyp = elem->point(n),
2065  nxpym = elem->point(n),
2066  nxpyp = elem->point(n);
2067  nxmym(0) -= TOLERANCE;
2068  nxmym(1) -= TOLERANCE;
2069  nxmyp(0) -= TOLERANCE;
2070  nxmyp(1) += TOLERANCE;
2071  nxpym(0) += TOLERANCE;
2072  nxpym(1) -= TOLERANCE;
2073  nxpyp(0) += TOLERANCE;
2074  nxpyp(1) += TOLERANCE;
2075  Gradient gxmym = g->component(context,
2076  var_component,
2077  nxmym,
2078  system.time);
2079  Gradient gxmyp = g->component(context,
2080  var_component,
2081  nxmyp,
2082  system.time);
2083  Gradient gxpym = g->component(context,
2084  var_component,
2085  nxpym,
2086  system.time);
2087  Gradient gxpyp = g->component(context,
2088  var_component,
2089  nxpyp,
2090  system.time);
2091  Number gxzplus = (gxpyp(2) - gxmyp(2))
2092  / 2. / TOLERANCE;
2093  Number gxzminus = (gxpym(2) - gxmym(2))
2094  / 2. / TOLERANCE;
2095  // xyz derivative
2096  Ue(current_dof) = (gxzplus - gxzminus)
2097  / 2. / TOLERANCE;
2098  dof_is_fixed[current_dof] = true;
2099  current_dof++;
2100  }
2101  }
2102  }
2103  // Assume that other C_ONE elements have a single nodal
2104  // value shape function and nodal gradient component
2105  // shape functions
2106  else if (cont == C_ONE)
2107  {
2108  libmesh_assert_equal_to (nc, 1 + dim);
2109  Ue(current_dof) = f->component(context,
2110  var_component,
2111  elem->point(n),
2112  system.time);
2113  dof_is_fixed[current_dof] = true;
2114  current_dof++;
2115  Gradient grad = g->component(context,
2116  var_component,
2117  elem->point(n),
2118  system.time);
2119  for (unsigned int i=0; i!= dim; ++i)
2120  {
2121  Ue(current_dof) = grad(i);
2122  dof_is_fixed[current_dof] = true;
2123  current_dof++;
2124  }
2125  }
2126  else
2127  libmesh_error();
2128  }
2129 
2130  // In 3D, project any edge values next
2131  if (dim > 2 && cont != DISCONTINUOUS)
2132  {
2133  const std::vector<Point>& xyz_values = edge_fe->get_xyz();
2134  const std::vector<Real>& JxW = edge_fe->get_JxW();
2135 
2136  const std::vector<std::vector<Real> >& phi = edge_fe->get_phi();
2137  const std::vector<std::vector<RealGradient> >* dphi = NULL;
2138  if (cont == C_ONE)
2139  dphi = &(edge_fe->get_dphi());
2140 
2141  for (unsigned int e=0; e != elem->n_edges(); ++e)
2142  {
2143  context.edge = e;
2144  context.edge_fe_reinit();
2145 
2146  const QBase& qedgerule = context.get_edge_qrule();
2147  const unsigned int n_qp = qedgerule.n_points();
2148 
2149  FEInterface::dofs_on_edge(elem, dim, fe_type, e,
2150  side_dofs);
2151 
2152  // Some edge dofs are on nodes and already
2153  // fixed, others are free to calculate
2154  unsigned int free_dofs = 0;
2155  for (unsigned int i=0; i != side_dofs.size(); ++i)
2156  if (!dof_is_fixed[side_dofs[i]])
2157  free_dof[free_dofs++] = i;
2158 
2159  // There may be nothing to project
2160  if (!free_dofs)
2161  continue;
2162 
2163  Ke.resize (free_dofs, free_dofs); Ke.zero();
2164  Fe.resize (free_dofs); Fe.zero();
2165  // The new edge coefficients
2166  DenseVector<Number> Uedge(free_dofs);
2167 
2168  // Loop over the quadrature points
2169  for (unsigned int qp=0; qp<n_qp; qp++)
2170  {
2171  // solution at the quadrature point
2172  Number fineval = f->component(context,
2173  var_component,
2174  xyz_values[qp],
2175  system.time);
2176  // solution grad at the quadrature point
2177  Gradient finegrad;
2178  if (cont == C_ONE)
2179  finegrad = g->component(context,
2180  var_component,
2181  xyz_values[qp],
2182  system.time);
2183 
2184  // Form edge projection matrix
2185  for (unsigned int sidei=0, freei=0;
2186  sidei != side_dofs.size(); ++sidei)
2187  {
2188  unsigned int i = side_dofs[sidei];
2189  // fixed DoFs aren't test functions
2190  if (dof_is_fixed[i])
2191  continue;
2192  for (unsigned int sidej=0, freej=0;
2193  sidej != side_dofs.size(); ++sidej)
2194  {
2195  unsigned int j = side_dofs[sidej];
2196  if (dof_is_fixed[j])
2197  Fe(freei) -= phi[i][qp] * phi[j][qp] *
2198  JxW[qp] * Ue(j);
2199  else
2200  Ke(freei,freej) += phi[i][qp] *
2201  phi[j][qp] * JxW[qp];
2202  if (cont == C_ONE)
2203  {
2204  if (dof_is_fixed[j])
2205  Fe(freei) -= ( (*dphi)[i][qp] *
2206  (*dphi)[j][qp] ) *
2207  JxW[qp] * Ue(j);
2208  else
2209  Ke(freei,freej) += ( (*dphi)[i][qp] *
2210  (*dphi)[j][qp] )
2211  * JxW[qp];
2212  }
2213  if (!dof_is_fixed[j])
2214  freej++;
2215  }
2216  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
2217  if (cont == C_ONE)
2218  Fe(freei) += (finegrad * (*dphi)[i][qp] ) *
2219  JxW[qp];
2220  freei++;
2221  }
2222  }
2223 
2224  Ke.cholesky_solve(Fe, Uedge);
2225 
2226  // Transfer new edge solutions to element
2227  for (unsigned int i=0; i != free_dofs; ++i)
2228  {
2229  Number &ui = Ue(side_dofs[free_dof[i]]);
2231  std::abs(ui - Uedge(i)) < TOLERANCE);
2232  ui = Uedge(i);
2233  dof_is_fixed[side_dofs[free_dof[i]]] = true;
2234  }
2235  }
2236  } // end if (dim > 2 && cont != DISCONTINUOUS)
2237 
2238  // Project any side values (edges in 2D, faces in 3D)
2239  if (dim > 1 && cont != DISCONTINUOUS)
2240  {
2241  const std::vector<Point>& xyz_values = side_fe->get_xyz();
2242  const std::vector<Real>& JxW = side_fe->get_JxW();
2243 
2244  const std::vector<std::vector<Real> >& phi = side_fe->get_phi();
2245  const std::vector<std::vector<RealGradient> >* dphi = NULL;
2246  if (cont == C_ONE)
2247  dphi = &(side_fe->get_dphi());
2248 
2249  for (unsigned int s=0; s != elem->n_sides(); ++s)
2250  {
2251  FEInterface::dofs_on_side(elem, dim, fe_type, s,
2252  side_dofs);
2253 
2254  // Some side dofs are on nodes/edges and already
2255  // fixed, others are free to calculate
2256  unsigned int free_dofs = 0;
2257  for (unsigned int i=0; i != side_dofs.size(); ++i)
2258  if (!dof_is_fixed[side_dofs[i]])
2259  free_dof[free_dofs++] = i;
2260 
2261  // There may be nothing to project
2262  if (!free_dofs)
2263  continue;
2264 
2265  Ke.resize (free_dofs, free_dofs); Ke.zero();
2266  Fe.resize (free_dofs); Fe.zero();
2267  // The new side coefficients
2268  DenseVector<Number> Uside(free_dofs);
2269 
2270  context.side = s;
2271  context.side_fe_reinit();
2272 
2273  const QBase& qsiderule = context.get_side_qrule();
2274  const unsigned int n_qp = qsiderule.n_points();
2275 
2276  // Loop over the quadrature points
2277  for (unsigned int qp=0; qp<n_qp; qp++)
2278  {
2279  // solution at the quadrature point
2280  Number fineval = f->component(context,
2281  var_component,
2282  xyz_values[qp],
2283  system.time);
2284  // solution grad at the quadrature point
2285  Gradient finegrad;
2286  if (cont == C_ONE)
2287  finegrad = g->component(context,
2288  var_component,
2289  xyz_values[qp],
2290  system.time);
2291 
2292  // Form side projection matrix
2293  for (unsigned int sidei=0, freei=0;
2294  sidei != side_dofs.size(); ++sidei)
2295  {
2296  unsigned int i = side_dofs[sidei];
2297  // fixed DoFs aren't test functions
2298  if (dof_is_fixed[i])
2299  continue;
2300  for (unsigned int sidej=0, freej=0;
2301  sidej != side_dofs.size(); ++sidej)
2302  {
2303  unsigned int j = side_dofs[sidej];
2304  if (dof_is_fixed[j])
2305  Fe(freei) -= phi[i][qp] * phi[j][qp] *
2306  JxW[qp] * Ue(j);
2307  else
2308  Ke(freei,freej) += phi[i][qp] *
2309  phi[j][qp] * JxW[qp];
2310  if (cont == C_ONE)
2311  {
2312  if (dof_is_fixed[j])
2313  Fe(freei) -= ( (*dphi)[i][qp] *
2314  (*dphi)[j][qp] ) *
2315  JxW[qp] * Ue(j);
2316  else
2317  Ke(freei,freej) += ( (*dphi)[i][qp] *
2318  (*dphi)[j][qp] )
2319  * JxW[qp];
2320  }
2321  if (!dof_is_fixed[j])
2322  freej++;
2323  }
2324  Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
2325  if (cont == C_ONE)
2326  Fe(freei) += (finegrad * (*dphi)[i][qp] ) *
2327  JxW[qp];
2328  freei++;
2329  }
2330  }
2331 
2332  Ke.cholesky_solve(Fe, Uside);
2333 
2334  // Transfer new side solutions to element
2335  for (unsigned int i=0; i != free_dofs; ++i)
2336  {
2337  Number &ui = Ue(side_dofs[free_dof[i]]);
2339  std::abs(ui - Uside(i)) < TOLERANCE);
2340  ui = Uside(i);
2341  dof_is_fixed[side_dofs[free_dof[i]]] = true;
2342  }
2343  }
2344  }// end if (dim > 1 && cont != DISCONTINUOUS)
2345 
2346  // Project the interior values, finally
2347 
2348  // Some interior dofs are on nodes/edges/sides and
2349  // already fixed, others are free to calculate
2350  unsigned int free_dofs = 0;
2351  for (unsigned int i=0; i != n_dofs; ++i)
2352  if (!dof_is_fixed[i])
2353  free_dof[free_dofs++] = i;
2354 
2355  // There may be nothing to project
2356  if (free_dofs)
2357  {
2358  context.elem_fe_reinit();
2359 
2360  const QBase& qrule = context.get_element_qrule();
2361  const unsigned int n_qp = qrule.n_points();
2362  const std::vector<Point>& xyz_values = fe->get_xyz();
2363  const std::vector<Real>& JxW = fe->get_JxW();
2364 
2365  const std::vector<std::vector<Real> >& phi = fe->get_phi();
2366  const std::vector<std::vector<RealGradient> >* dphi = NULL;
2367  if (cont == C_ONE)
2368  dphi = &(side_fe->get_dphi());
2369 
2370  Ke.resize (free_dofs, free_dofs); Ke.zero();
2371  Fe.resize (free_dofs); Fe.zero();
2372  // The new interior coefficients
2373  DenseVector<Number> Uint(free_dofs);
2374 
2375  // Loop over the quadrature points
2376  for (unsigned int qp=0; qp<n_qp; qp++)
2377  {
2378  // solution at the quadrature point
2379  Number fineval = f->component(context,
2380  var_component,
2381  xyz_values[qp],
2382  system.time);
2383  // solution grad at the quadrature point
2384  Gradient finegrad;
2385  if (cont == C_ONE)
2386  finegrad = g->component(context,
2387  var_component,
2388  xyz_values[qp],
2389  system.time);
2390 
2391  // Form interior projection matrix
2392  for (unsigned int i=0, freei=0; i != n_dofs; ++i)
2393  {
2394  // fixed DoFs aren't test functions
2395  if (dof_is_fixed[i])
2396  continue;
2397  for (unsigned int j=0, freej=0; j != n_dofs; ++j)
2398  {
2399  if (dof_is_fixed[j])
2400  Fe(freei) -= phi[i][qp] * phi[j][qp] * JxW[qp]
2401  * Ue(j);
2402  else
2403  Ke(freei,freej) += phi[i][qp] * phi[j][qp] *
2404  JxW[qp];
2405  if (cont == C_ONE)
2406  {
2407  if (dof_is_fixed[j])
2408  Fe(freei) -= ( (*dphi)[i][qp] *
2409  (*dphi)[j][qp] ) * JxW[qp] *
2410  Ue(j);
2411  else
2412  Ke(freei,freej) += ( (*dphi)[i][qp] *
2413  (*dphi)[j][qp] ) *
2414  JxW[qp];
2415  }
2416  if (!dof_is_fixed[j])
2417  freej++;
2418  }
2419  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
2420  if (cont == C_ONE)
2421  Fe(freei) += (finegrad * (*dphi)[i][qp] ) * JxW[qp];
2422  freei++;
2423  }
2424  }
2425  Ke.cholesky_solve(Fe, Uint);
2426 
2427  // Transfer new interior solutions to element
2428  for (unsigned int i=0; i != free_dofs; ++i)
2429  {
2430  Number &ui = Ue(free_dof[i]);
2432  std::abs(ui - Uint(i)) < TOLERANCE);
2433  ui = Uint(i);
2434  dof_is_fixed[free_dof[i]] = true;
2435  }
2436 
2437  } // if there are free interior dofs
2438 
2439  // Make sure every DoF got reached!
2440  for (unsigned int i=0; i != n_dofs; ++i)
2441  libmesh_assert(dof_is_fixed[i]);
2442 
2443  const numeric_index_type
2444  first = new_vector.first_local_index(),
2445  last = new_vector.last_local_index();
2446 
2447  // Lock the new_vector since it is shared among threads.
2448  {
2449  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
2450 
2451  for (unsigned int i = 0; i < n_dofs; i++)
2452  // We may be projecting a new zero value onto
2453  // an old nonzero approximation - RHS
2454  // if (Ue(i) != 0.)
2455  if ((dof_indices[i] >= first) &&
2456  (dof_indices[i] < last))
2457  {
2458  new_vector.set(dof_indices[i], Ue(i));
2459  }
2460  }
2461  } // end variables loop
2462  } // end elem loop
2463 }

Member Data Documentation

AutoPtr<FEMFunctionBase<Number> > libMesh::ProjectFEMSolution::f
private

Definition at line 165 of file system_projection.C.

Referenced by ProjectFEMSolution().

AutoPtr<FEMFunctionBase<Gradient> > libMesh::ProjectFEMSolution::g
private

Definition at line 166 of file system_projection.C.

NumericVector<Number>& libMesh::ProjectFEMSolution::new_vector
private

Definition at line 167 of file system_projection.C.

const System& libMesh::ProjectFEMSolution::system
private

Definition at line 163 of file system_projection.C.

Referenced by operator()().


The documentation for this class was generated from the following file:

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:08 UTC

Hosted By:
SourceForge.net Logo