libMesh::PatchRecoveryErrorEstimator::EstimateError Class Reference
Public Member Functions | |
| EstimateError (const System &sys, const PatchRecoveryErrorEstimator &ee, ErrorVector &epc) | |
| void | operator() (const ConstElemRange &range) const |
Private Attributes | |
| const System & | system |
| const PatchRecoveryErrorEstimator & | error_estimator |
| ErrorVector & | error_per_cell |
Detailed Description
Class to compute the error contribution for a range of elements. May be executed in parallel on separate threads.
Definition at line 114 of file patch_recovery_error_estimator.h.
Constructor & Destructor Documentation
| libMesh::PatchRecoveryErrorEstimator::EstimateError::EstimateError | ( | const System & | sys, | |
| const PatchRecoveryErrorEstimator & | ee, | |||
| ErrorVector & | epc | |||
| ) | [inline] |
Definition at line 117 of file patch_recovery_error_estimator.h.
00119 : 00120 system(sys), 00121 error_estimator(ee), 00122 error_per_cell(epc) 00123 {}
Member Function Documentation
| void libMesh::PatchRecoveryErrorEstimator::EstimateError::operator() | ( | const ConstElemRange & | range | ) | const |
Definition at line 188 of file patch_recovery_error_estimator.C.
References std::abs(), libMesh::TypeTensor< T >::add_scaled(), libMesh::TypeVector< T >::add_scaled(), libMesh::StoredRange< iterator_type, object_type >::begin(), libMesh::FEGenericBase< Real >::build(), libMesh::Patch::build_around_element(), libMesh::System::current_solution(), libMesh::FEType::default_quadrature_rule(), libMesh::DofMap::dof_indices(), libMesh::StoredRange< iterator_type, object_type >::end(), libMesh::err, error_estimator, libMesh::ErrorEstimator::error_norm, error_per_cell, libMesh::ErrorVectorReal, libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMeshEnums::H1_SEMINORM, libMeshEnums::H1_X_SEMINORM, libMeshEnums::H1_Y_SEMINORM, libMeshEnums::H1_Z_SEMINORM, libMeshEnums::H2_SEMINORM, libMesh::DofObject::id(), libMeshEnums::L2, libMeshEnums::L_INF, libMesh::DenseMatrix< T >::lu_solve(), libMesh::DenseMatrixBase< T >::m(), std::max(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::DenseMatrixBase< T >::n(), libMesh::System::n_vars(), n_vars, libMesh::TensorTools::norm_sq(), libMesh::FEType::order, libMesh::Elem::p_level(), libMesh::PatchRecoveryErrorEstimator::patch_growth_strategy, libMesh::PatchRecoveryErrorEstimator::patch_reuse, libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::PatchRecoveryErrorEstimator::specpoly(), libMesh::Threads::spin_mtx, system, libMesh::PatchRecoveryErrorEstimator::target_patch_size, libMesh::SystemNorm::type(), libMesh::DofMap::variable_type(), libMeshEnums::W1_INF_SEMINORM, libMeshEnums::W2_INF_SEMINORM, libMesh::SystemNorm::weight(), libMesh::SystemNorm::weight_sq(), and libMesh::zero.
00189 { 00190 // The current mesh 00191 const MeshBase& mesh = system.get_mesh(); 00192 00193 // The dimensionality of the mesh 00194 const unsigned int dim = mesh.mesh_dimension(); 00195 00196 // The number of variables in the system 00197 const unsigned int n_vars = system.n_vars(); 00198 00199 // The DofMap for this system 00200 const DofMap& dof_map = system.get_dof_map(); 00201 00202 //------------------------------------------------------------ 00203 // Iterate over all the elements in the range. 00204 for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it!=range.end(); ++elem_it) 00205 { 00206 // elem is necessarily an active element on the local processor 00207 const Elem* elem = *elem_it; 00208 00209 // We'll need an index into the error vector 00210 const dof_id_type e_id=elem->id(); 00211 00212 // We are going to build a patch containing the current element 00213 // and its neighbors on the local processor 00214 Patch patch; 00215 00216 // If we are reusing patches and the current element 00217 // already has an estimate associated with it, move on the 00218 // next element 00219 if(this->error_estimator.patch_reuse && error_per_cell[e_id] != 0) 00220 continue; 00221 00222 // If we are not reusing patches or havent built one containing this element, we build one 00223 00224 // Use user specified patch size and growth strategy 00225 patch.build_around_element (elem, error_estimator.target_patch_size, 00226 error_estimator.patch_growth_strategy); 00227 00228 // Declare a new_error_per_cell vector to hold error estimates 00229 // from each element in this patch, or one estimate if we are 00230 // not reusing patches since we will only be computing error for 00231 // one cell 00232 std::vector<Real> new_error_per_cell(1, 0.); 00233 if(this->error_estimator.patch_reuse) 00234 new_error_per_cell.resize(patch.size(), 0.); 00235 00236 //------------------------------------------------------------ 00237 // Process each variable in the system using the current patch 00238 for (unsigned int var=0; var<n_vars; var++) 00239 { 00240 #ifndef DEBUG 00241 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 00242 libmesh_assert (error_estimator.error_norm.type(var) == L2 || 00243 error_estimator.error_norm.type(var) == H1_SEMINORM || 00244 error_estimator.error_norm.type(var) == H2_SEMINORM || 00245 error_estimator.error_norm.type(var) == H1_X_SEMINORM || 00246 error_estimator.error_norm.type(var) == H1_Y_SEMINORM || 00247 error_estimator.error_norm.type(var) == H1_Z_SEMINORM || 00248 error_estimator.error_norm.type(var) == L_INF || 00249 error_estimator.error_norm.type(var) == W1_INF_SEMINORM || 00250 error_estimator.error_norm.type(var) == W2_INF_SEMINORM); 00251 #else 00252 libmesh_assert (error_estimator.error_norm.type(var) == L2 || 00253 error_estimator.error_norm.type(var) == L_INF || 00254 error_estimator.error_norm.type(var) == H1_SEMINORM || 00255 error_estimator.error_norm.type(var) == H1_X_SEMINORM || 00256 error_estimator.error_norm.type(var) == H1_Y_SEMINORM || 00257 error_estimator.error_norm.type(var) == H1_Z_SEMINORM || 00258 error_estimator.error_norm.type(var) == W1_INF_SEMINORM); 00259 #endif 00260 if (var > 0) 00261 // We can't mix L_inf and L_2 norms 00262 libmesh_assert (((error_estimator.error_norm.type(var) == L2 || 00263 error_estimator.error_norm.type(var) == H1_SEMINORM || 00264 error_estimator.error_norm.type(var) == H1_X_SEMINORM || 00265 error_estimator.error_norm.type(var) == H1_Y_SEMINORM || 00266 error_estimator.error_norm.type(var) == H1_Z_SEMINORM || 00267 error_estimator.error_norm.type(var) == H2_SEMINORM) && 00268 (error_estimator.error_norm.type(var-1) == L2 || 00269 error_estimator.error_norm.type(var-1) == H1_SEMINORM || 00270 error_estimator.error_norm.type(var-1) == H1_X_SEMINORM || 00271 error_estimator.error_norm.type(var-1) == H1_Y_SEMINORM || 00272 error_estimator.error_norm.type(var-1) == H1_Z_SEMINORM || 00273 error_estimator.error_norm.type(var-1) == H2_SEMINORM)) || 00274 ((error_estimator.error_norm.type(var) == L_INF || 00275 error_estimator.error_norm.type(var) == W1_INF_SEMINORM || 00276 error_estimator.error_norm.type(var) == W2_INF_SEMINORM) && 00277 (error_estimator.error_norm.type(var-1) == L_INF || 00278 error_estimator.error_norm.type(var-1) == W1_INF_SEMINORM || 00279 error_estimator.error_norm.type(var-1) == W2_INF_SEMINORM))); 00280 #endif 00281 00282 // Possibly skip this variable 00283 if (error_estimator.error_norm.weight(var) == 0.0) continue; 00284 00285 // The type of finite element to use for this variable 00286 const FEType& fe_type = dof_map.variable_type (var); 00287 00288 const Order element_order = static_cast<Order> 00289 (fe_type.order + elem->p_level()); 00290 00291 // Finite element object for use in this patch 00292 AutoPtr<FEBase> fe (FEBase::build (dim, fe_type)); 00293 00294 // Build an appropriate Gaussian quadrature rule 00295 AutoPtr<QBase> qrule (fe_type.default_quadrature_rule(dim)); 00296 00297 // Tell the finite element about the quadrature rule. 00298 fe->attach_quadrature_rule (qrule.get()); 00299 00300 // Get Jacobian values, etc.. 00301 const std::vector<Real>& JxW = fe->get_JxW(); 00302 const std::vector<Point>& q_point = fe->get_xyz(); 00303 00304 // Get whatever phi/dphi/d2phi values we need. Avoid 00305 // getting them unless the requested norm is actually going 00306 // to use them. 00307 00308 const std::vector<std::vector<Real> > *phi = NULL; 00309 // If we're using phi to assert the correct dof_indices 00310 // vector size later, then we'll need to get_phi whether we 00311 // plan to use it or not. 00312 #ifdef NDEBUG 00313 if (error_estimator.error_norm.type(var) == L2 || 00314 error_estimator.error_norm.type(var) == L_INF) 00315 #endif 00316 phi = &(fe->get_phi()); 00317 00318 const std::vector<std::vector<RealGradient> > *dphi = NULL; 00319 if (error_estimator.error_norm.type(var) == H1_SEMINORM || 00320 error_estimator.error_norm.type(var) == H1_X_SEMINORM || 00321 error_estimator.error_norm.type(var) == H1_Y_SEMINORM || 00322 error_estimator.error_norm.type(var) == H1_Z_SEMINORM || 00323 error_estimator.error_norm.type(var) == W1_INF_SEMINORM) 00324 dphi = &(fe->get_dphi()); 00325 00326 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 00327 const std::vector<std::vector<RealTensor> > *d2phi = NULL; 00328 if (error_estimator.error_norm.type(var) == H2_SEMINORM || 00329 error_estimator.error_norm.type(var) == W2_INF_SEMINORM) 00330 d2phi = &(fe->get_d2phi()); 00331 #endif 00332 00333 // global DOF indices 00334 std::vector<dof_id_type> dof_indices; 00335 00336 // Compute the approprite size for the patch projection matrices 00337 // and vectors; 00338 unsigned int matsize = element_order + 1; 00339 if (dim > 1) 00340 { 00341 matsize *= (element_order + 2); 00342 matsize /= 2; 00343 } 00344 if (dim > 2) 00345 { 00346 matsize *= (element_order + 3); 00347 matsize /= 3; 00348 } 00349 00350 DenseMatrix<Number> Kp(matsize,matsize); 00351 DenseVector<Number> F, Fx, Fy, Fz, Fxy, Fxz, Fyz; 00352 DenseVector<Number> Pu_h, Pu_x_h, Pu_y_h, Pu_z_h, Pu_xy_h, Pu_xz_h, Pu_yz_h; 00353 if (error_estimator.error_norm.type(var) == L2 || 00354 error_estimator.error_norm.type(var) == L_INF) 00355 { 00356 F.resize(matsize); Pu_h.resize(matsize); 00357 } 00358 else if (error_estimator.error_norm.type(var) == H1_SEMINORM || 00359 error_estimator.error_norm.type(var) == W1_INF_SEMINORM || 00360 error_estimator.error_norm.type(var) == H2_SEMINORM || 00361 error_estimator.error_norm.type(var) == W2_INF_SEMINORM) 00362 { 00363 Fx.resize(matsize); Pu_x_h.resize(matsize); // stores xx in W2 cases 00364 #if LIBMESH_DIM > 1 00365 Fy.resize(matsize); Pu_y_h.resize(matsize); // stores yy in W2 cases 00366 #endif 00367 #if LIBMESH_DIM > 2 00368 Fz.resize(matsize); Pu_z_h.resize(matsize); // stores zz in W2 cases 00369 #endif 00370 } 00371 else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM) 00372 { 00373 Fx.resize(matsize); Pu_x_h.resize(matsize); // Only need to compute the x gradient for the x component seminorm 00374 } 00375 else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM) 00376 { 00377 libmesh_assert_greater (LIBMESH_DIM, 1); 00378 Fy.resize(matsize); Pu_y_h.resize(matsize); // Only need to compute the y gradient for the y component seminorm 00379 } 00380 else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM) 00381 { 00382 libmesh_assert_greater (LIBMESH_DIM, 2); 00383 Fz.resize(matsize); Pu_z_h.resize(matsize); // Only need to compute the z gradient for the z component seminorm 00384 } 00385 00386 #if LIBMESH_DIM > 1 00387 if (error_estimator.error_norm.type(var) == H2_SEMINORM || 00388 error_estimator.error_norm.type(var) == W2_INF_SEMINORM) 00389 { 00390 Fxy.resize(matsize); Pu_xy_h.resize(matsize); 00391 #if LIBMESH_DIM > 2 00392 Fxz.resize(matsize); Pu_xz_h.resize(matsize); 00393 Fyz.resize(matsize); Pu_yz_h.resize(matsize); 00394 #endif 00395 } 00396 #endif 00397 00398 //------------------------------------------------------ 00399 // Loop over each element in the patch and compute their 00400 // contribution to the patch gradient projection. 00401 Patch::const_iterator patch_it = patch.begin(); 00402 const Patch::const_iterator patch_end = patch.end(); 00403 00404 for (; patch_it != patch_end; ++patch_it) 00405 { 00406 // The pth element in the patch 00407 const Elem* e_p = *patch_it; 00408 00409 // Reinitialize the finite element data for this element 00410 fe->reinit (e_p); 00411 00412 // Get the global DOF indices for the current variable 00413 // in the current element 00414 dof_map.dof_indices (e_p, dof_indices, var); 00415 libmesh_assert_equal_to (dof_indices.size(), phi->size()); 00416 00417 const unsigned int n_dofs = 00418 libmesh_cast_int<unsigned int>(dof_indices.size()); 00419 const unsigned int n_qp = qrule->n_points(); 00420 00421 // Compute the projection components from this cell. 00422 // \int_{Omega_e} \psi_i \psi_j = \int_{Omega_e} du_h/dx_k \psi_i 00423 for (unsigned int qp=0; qp<n_qp; qp++) 00424 { 00425 // Construct the shape function values for the patch projection 00426 std::vector<Real> psi(specpoly(dim, element_order, q_point[qp], matsize)); 00427 00428 // Patch matrix contribution 00429 for (unsigned int i=0; i<Kp.m(); i++) 00430 for (unsigned int j=0; j<Kp.n(); j++) 00431 Kp(i,j) += JxW[qp]*psi[i]*psi[j]; 00432 00433 if (error_estimator.error_norm.type(var) == L2 || 00434 error_estimator.error_norm.type(var) == L_INF) 00435 { 00436 // Compute the solution on the current patch element 00437 // the quadrature point 00438 Number u_h = libMesh::zero; 00439 00440 for (unsigned int i=0; i<n_dofs; i++) 00441 u_h += (*phi)[i][qp]*system.current_solution (dof_indices[i]); 00442 00443 // Patch RHS contributions 00444 for (unsigned int i=0; i<psi.size(); i++) 00445 F(i) = JxW[qp]*u_h*psi[i]; 00446 00447 } 00448 else if (error_estimator.error_norm.type(var) == H1_SEMINORM || 00449 error_estimator.error_norm.type(var) == W1_INF_SEMINORM) 00450 { 00451 // Compute the gradient on the current patch element 00452 // at the quadrature point 00453 Gradient grad_u_h; 00454 00455 for (unsigned int i=0; i<n_dofs; i++) 00456 grad_u_h.add_scaled ((*dphi)[i][qp], 00457 system.current_solution(dof_indices[i])); 00458 00459 // Patch RHS contributions 00460 for (unsigned int i=0; i<psi.size(); i++) 00461 { 00462 Fx(i) += JxW[qp]*grad_u_h(0)*psi[i]; 00463 #if LIBMESH_DIM > 1 00464 Fy(i) += JxW[qp]*grad_u_h(1)*psi[i]; 00465 #endif 00466 #if LIBMESH_DIM > 2 00467 Fz(i) += JxW[qp]*grad_u_h(2)*psi[i]; 00468 #endif 00469 } 00470 } 00471 else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM) 00472 { 00473 // Compute the gradient on the current patch element 00474 // at the quadrature point 00475 Gradient grad_u_h; 00476 00477 for (unsigned int i=0; i<n_dofs; i++) 00478 grad_u_h.add_scaled ((*dphi)[i][qp], 00479 system.current_solution(dof_indices[i])); 00480 00481 // Patch RHS contributions 00482 for (unsigned int i=0; i<psi.size(); i++) 00483 { 00484 Fx(i) += JxW[qp]*grad_u_h(0)*psi[i]; 00485 } 00486 } 00487 else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM) 00488 { 00489 // Compute the gradient on the current patch element 00490 // at the quadrature point 00491 Gradient grad_u_h; 00492 00493 for (unsigned int i=0; i<n_dofs; i++) 00494 grad_u_h.add_scaled ((*dphi)[i][qp], 00495 system.current_solution(dof_indices[i])); 00496 00497 // Patch RHS contributions 00498 for (unsigned int i=0; i<psi.size(); i++) 00499 { 00500 Fy(i) += JxW[qp]*grad_u_h(1)*psi[i]; 00501 } 00502 } 00503 else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM) 00504 { 00505 // Compute the gradient on the current patch element 00506 // at the quadrature point 00507 Gradient grad_u_h; 00508 00509 for (unsigned int i=0; i<n_dofs; i++) 00510 grad_u_h.add_scaled ((*dphi)[i][qp], 00511 system.current_solution(dof_indices[i])); 00512 00513 // Patch RHS contributions 00514 for (unsigned int i=0; i<psi.size(); i++) 00515 { 00516 Fz(i) += JxW[qp]*grad_u_h(2)*psi[i]; 00517 } 00518 } 00519 else if (error_estimator.error_norm.type(var) == H2_SEMINORM || 00520 error_estimator.error_norm.type(var) == W2_INF_SEMINORM) 00521 { 00522 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 00523 // Compute the hessian on the current patch element 00524 // at the quadrature point 00525 Tensor hess_u_h; 00526 00527 for (unsigned int i=0; i<n_dofs; i++) 00528 hess_u_h.add_scaled ((*d2phi)[i][qp], 00529 system.current_solution(dof_indices[i])); 00530 00531 // Patch RHS contributions 00532 for (unsigned int i=0; i<psi.size(); i++) 00533 { 00534 Fx(i) += JxW[qp]*hess_u_h(0,0)*psi[i]; 00535 #if LIBMESH_DIM > 1 00536 Fy(i) += JxW[qp]*hess_u_h(1,1)*psi[i]; 00537 Fxy(i) += JxW[qp]*hess_u_h(0,1)*psi[i]; 00538 #endif 00539 #if LIBMESH_DIM > 2 00540 Fz(i) += JxW[qp]*hess_u_h(2,2)*psi[i]; 00541 Fxz(i) += JxW[qp]*hess_u_h(0,2)*psi[i]; 00542 Fyz(i) += JxW[qp]*hess_u_h(1,2)*psi[i]; 00543 #endif 00544 } 00545 #else 00546 libMesh::err << "ERROR: --enable-second-derivatives is required\n" 00547 << " for _sobolev_order == 2!\n"; 00548 libmesh_error(); 00549 #endif 00550 } 00551 else 00552 libmesh_error(); 00553 } // end quadrature loop 00554 } // end patch loop 00555 00556 00557 00558 //-------------------------------------------------- 00559 // Now we have fully assembled the projection system 00560 // for this patch. Project the gradient components. 00561 // MAY NEED TO USE PARTIAL PIVOTING! 00562 if (error_estimator.error_norm.type(var) == L2 || 00563 error_estimator.error_norm.type(var) == L_INF) 00564 { 00565 Kp.lu_solve(F, Pu_h); 00566 } 00567 else if (error_estimator.error_norm.type(var) == H1_SEMINORM || 00568 error_estimator.error_norm.type(var) == W1_INF_SEMINORM || 00569 error_estimator.error_norm.type(var) == H2_SEMINORM || 00570 error_estimator.error_norm.type(var) == W2_INF_SEMINORM) 00571 { 00572 Kp.lu_solve (Fx, Pu_x_h); 00573 #if LIBMESH_DIM > 1 00574 Kp.lu_solve (Fy, Pu_y_h); 00575 #endif 00576 #if LIBMESH_DIM > 2 00577 Kp.lu_solve (Fz, Pu_z_h); 00578 #endif 00579 } 00580 else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM) 00581 { 00582 Kp.lu_solve (Fx, Pu_x_h); 00583 } 00584 else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM) 00585 { 00586 Kp.lu_solve (Fy, Pu_y_h); 00587 } 00588 else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM) 00589 { 00590 Kp.lu_solve (Fz, Pu_z_h); 00591 } 00592 00593 #if LIBMESH_DIM > 1 00594 if (error_estimator.error_norm.type(var) == H2_SEMINORM || 00595 error_estimator.error_norm.type(var) == W2_INF_SEMINORM) 00596 { 00597 Kp.lu_solve(Fxy, Pu_xy_h); 00598 #if LIBMESH_DIM > 2 00599 Kp.lu_solve(Fxz, Pu_xz_h); 00600 Kp.lu_solve(Fyz, Pu_yz_h); 00601 #endif 00602 } 00603 #endif 00604 00605 // If we are reusing patches, reuse the current patch to loop 00606 // over all elements in the current patch, otherwise build a new 00607 // patch containing just the current element and loop over it 00608 // Note that C++ will not allow patch_re_end to be a const here 00609 Patch::const_iterator patch_re_it; 00610 Patch::const_iterator patch_re_end; 00611 00612 // Declare a new patch 00613 Patch patch_re; 00614 00615 if(this->error_estimator.patch_reuse) 00616 { 00617 // Just get the iterators from the current patch 00618 patch_re_it = patch.begin(); 00619 patch_re_end = patch.end(); 00620 } 00621 else 00622 { 00623 // Use a target patch size of just 0, this will contain 00624 // just the current element 00625 patch_re.build_around_element (elem, 0, 00626 error_estimator.patch_growth_strategy); 00627 00628 // Get the iterators from this newly constructed patch 00629 patch_re_it = patch_re.begin(); 00630 patch_re_end = patch_re.end(); 00631 } 00632 00633 // If we are reusing patches, loop over all the elements 00634 // in the current patch and develop an estimate 00635 // for all the elements by computing ||P u_h - u_h|| or ||P grad_u_h - grad_u_h|| 00636 // or ||P hess_u_h - hess_u_h|| according to the requested 00637 // seminorm, otherwise just compute it for the current element 00638 00639 // Loop over every element in the patch 00640 for (unsigned int e = 0 ; patch_re_it != patch_re_end; patch_re_it++, ++e) 00641 { 00642 // Build the Finite Element for the current element 00643 00644 // The pth element in the patch 00645 const Elem* e_p = *patch_re_it; 00646 00647 // We'll need an index into the error vector for this element 00648 const dof_id_type e_p_id = e_p->id(); 00649 00650 // We will update the new_error_per_cell vector with element_error if the 00651 // error_per_cell[e_p_id] entry is non-zero, otherwise update it 00652 // with 0. i.e. leave it unchanged 00653 00654 // No need to compute the estimate if we are reusing patches and already have one 00655 if (this->error_estimator.patch_reuse && error_per_cell[e_p_id] != 0.) 00656 continue; 00657 00658 // Reinitialize the finite element data for this element 00659 fe->reinit (e_p); 00660 00661 // Get the global DOF indices for the current variable 00662 // in the current element 00663 dof_map.dof_indices (e_p, dof_indices, var); 00664 libmesh_assert_equal_to (dof_indices.size(), phi->size()); 00665 00666 // The number of dofs for this variable on this element 00667 const unsigned int n_dofs = 00668 libmesh_cast_int<unsigned int>(dof_indices.size()); 00669 00670 // Variable to hold the error on the current element 00671 Real element_error = 0; 00672 00673 const Order qorder = 00674 static_cast<Order>(fe_type.order + e_p->p_level()); 00675 00676 // A quadrature rule for this element 00677 QGrid samprule (dim, qorder); 00678 00679 if (error_estimator.error_norm.type(var) == W1_INF_SEMINORM || 00680 error_estimator.error_norm.type(var) == W2_INF_SEMINORM) 00681 fe->attach_quadrature_rule (&samprule); 00682 00683 // The number of points we will sample over 00684 const unsigned int n_sp = 00685 libmesh_cast_int<unsigned int>(JxW.size()); 00686 00687 // Loop over every sample point for the current element 00688 for (unsigned int sp=0; sp<n_sp; sp++) 00689 { 00690 // Compute the solution at the current sample point 00691 00692 std::vector<Number> temperr(6,0.0); // x,y,z or xx,yy,zz,xy,xz,yz 00693 00694 if (error_estimator.error_norm.type(var) == L2 || 00695 error_estimator.error_norm.type(var) == L_INF) 00696 { 00697 // Compute the value at the current sample point 00698 Number u_h = libMesh::zero; 00699 00700 for (unsigned int i=0; i<n_dofs; i++) 00701 u_h += (*phi)[i][sp]*system.current_solution (dof_indices[i]); 00702 00703 // Compute the phi values at the current sample point 00704 std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize)); 00705 for (unsigned int i=0; i<matsize; i++) 00706 { 00707 temperr[0] += psi[i]*Pu_h(i); 00708 } 00709 00710 temperr[0] -= u_h; 00711 } 00712 else if (error_estimator.error_norm.type(var) == H1_SEMINORM || 00713 error_estimator.error_norm.type(var) == W1_INF_SEMINORM) 00714 { 00715 // Compute the gradient at the current sample point 00716 Gradient grad_u_h; 00717 00718 for (unsigned int i=0; i<n_dofs; i++) 00719 grad_u_h.add_scaled ((*dphi)[i][sp], 00720 system.current_solution(dof_indices[i])); 00721 00722 // Compute the phi values at the current sample point 00723 std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize)); 00724 00725 for (unsigned int i=0; i<matsize; i++) 00726 { 00727 temperr[0] += psi[i]*Pu_x_h(i); 00728 #if LIBMESH_DIM > 1 00729 temperr[1] += psi[i]*Pu_y_h(i); 00730 #endif 00731 #if LIBMESH_DIM > 2 00732 temperr[2] += psi[i]*Pu_z_h(i); 00733 #endif 00734 } 00735 temperr[0] -= grad_u_h(0); 00736 #if LIBMESH_DIM > 1 00737 temperr[1] -= grad_u_h(1); 00738 #endif 00739 #if LIBMESH_DIM > 2 00740 temperr[2] -= grad_u_h(2); 00741 #endif 00742 } 00743 else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM) 00744 { 00745 // Compute the gradient at the current sample point 00746 Gradient grad_u_h; 00747 00748 for (unsigned int i=0; i<n_dofs; i++) 00749 grad_u_h.add_scaled ((*dphi)[i][sp], 00750 system.current_solution(dof_indices[i])); 00751 00752 // Compute the phi values at the current sample point 00753 std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize)); 00754 for (unsigned int i=0; i<matsize; i++) 00755 { 00756 temperr[0] += psi[i]*Pu_x_h(i); 00757 } 00758 00759 temperr[0] -= grad_u_h(0); 00760 } 00761 else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM) 00762 { 00763 // Compute the gradient at the current sample point 00764 Gradient grad_u_h; 00765 00766 for (unsigned int i=0; i<n_dofs; i++) 00767 grad_u_h.add_scaled ((*dphi)[i][sp], 00768 system.current_solution(dof_indices[i])); 00769 00770 // Compute the phi values at the current sample point 00771 std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize)); 00772 for (unsigned int i=0; i<matsize; i++) 00773 { 00774 temperr[1] += psi[i]*Pu_y_h(i); 00775 } 00776 00777 temperr[1] -= grad_u_h(1); 00778 } 00779 else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM) 00780 { 00781 // Compute the gradient at the current sample point 00782 Gradient grad_u_h; 00783 00784 for (unsigned int i=0; i<n_dofs; i++) 00785 grad_u_h.add_scaled ((*dphi)[i][sp], 00786 system.current_solution(dof_indices[i])); 00787 00788 // Compute the phi values at the current sample point 00789 std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize)); 00790 for (unsigned int i=0; i<matsize; i++) 00791 { 00792 temperr[2] += psi[i]*Pu_z_h(i); 00793 } 00794 00795 temperr[2] -= grad_u_h(2); 00796 } 00797 else if (error_estimator.error_norm.type(var) == H2_SEMINORM || 00798 error_estimator.error_norm.type(var) == W2_INF_SEMINORM) 00799 { 00800 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 00801 // Compute the Hessian at the current sample point 00802 Tensor hess_u_h; 00803 00804 for (unsigned int i=0; i<n_dofs; i++) 00805 hess_u_h.add_scaled ((*d2phi)[i][sp], 00806 system.current_solution(dof_indices[i])); 00807 00808 // Compute the phi values at the current sample point 00809 std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize)); 00810 for (unsigned int i=0; i<matsize; i++) 00811 { 00812 temperr[0] += psi[i]*Pu_x_h(i); 00813 #if LIBMESH_DIM > 1 00814 temperr[1] += psi[i]*Pu_y_h(i); 00815 temperr[3] += psi[i]*Pu_xy_h(i); 00816 #endif 00817 #if LIBMESH_DIM > 2 00818 temperr[2] += psi[i]*Pu_z_h(i); 00819 temperr[4] += psi[i]*Pu_xz_h(i); 00820 temperr[5] += psi[i]*Pu_yz_h(i); 00821 #endif 00822 } 00823 00824 temperr[0] -= hess_u_h(0,0); 00825 #if LIBMESH_DIM > 1 00826 temperr[1] -= hess_u_h(1,1); 00827 temperr[3] -= hess_u_h(0,1); 00828 #endif 00829 #if LIBMESH_DIM > 2 00830 temperr[2] -= hess_u_h(2,2); 00831 temperr[4] -= hess_u_h(0,2); 00832 temperr[5] -= hess_u_h(1,2); 00833 #endif 00834 #else 00835 libMesh::err << "ERROR: --enable-second-derivatives is required\n" 00836 << " for _sobolev_order == 2!\n"; 00837 libmesh_error(); 00838 #endif 00839 } 00840 // Add up relevant terms. We can easily optimize the 00841 // LIBMESH_DIM < 3 cases a little bit with the exception 00842 // of the W2 cases 00843 00844 if (error_estimator.error_norm.type(var) == L_INF) 00845 element_error = std::max(element_error, std::abs(temperr[0])); 00846 else if (error_estimator.error_norm.type(var) == W1_INF_SEMINORM) 00847 for (unsigned int i=0; i != LIBMESH_DIM; ++i) 00848 element_error = std::max(element_error, std::abs(temperr[i])); 00849 else if (error_estimator.error_norm.type(var) == W2_INF_SEMINORM) 00850 for (unsigned int i=0; i != 6; ++i) 00851 element_error = std::max(element_error, std::abs(temperr[i])); 00852 else if (error_estimator.error_norm.type(var) == L2) 00853 element_error += JxW[sp]*TensorTools::norm_sq(temperr[0]); 00854 else if (error_estimator.error_norm.type(var) == H1_SEMINORM) 00855 for (unsigned int i=0; i != LIBMESH_DIM; ++i) 00856 element_error += JxW[sp]*TensorTools::norm_sq(temperr[i]); 00857 else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM) 00858 element_error += JxW[sp]*TensorTools::norm_sq(temperr[0]); 00859 else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM) 00860 element_error += JxW[sp]*TensorTools::norm_sq(temperr[1]); 00861 else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM) 00862 element_error += JxW[sp]*TensorTools::norm_sq(temperr[2]); 00863 else if (error_estimator.error_norm.type(var) == H2_SEMINORM) 00864 { 00865 for (unsigned int i=0; i != LIBMESH_DIM; ++i) 00866 element_error += JxW[sp]*TensorTools::norm_sq(temperr[i]); 00867 // Off diagonal terms enter into the Hessian norm twice 00868 for (unsigned int i=3; i != 6; ++i) 00869 element_error += JxW[sp]*2*TensorTools::norm_sq(temperr[i]); 00870 } 00871 00872 } // End loop over sample points 00873 00874 if (error_estimator.error_norm.type(var) == L_INF || 00875 error_estimator.error_norm.type(var) == W1_INF_SEMINORM || 00876 error_estimator.error_norm.type(var) == W2_INF_SEMINORM) 00877 new_error_per_cell[e] += error_estimator.error_norm.weight(var) * element_error; 00878 else if (error_estimator.error_norm.type(var) == L2 || 00879 error_estimator.error_norm.type(var) == H1_SEMINORM || 00880 error_estimator.error_norm.type(var) == H1_X_SEMINORM || 00881 error_estimator.error_norm.type(var) == H1_Y_SEMINORM || 00882 error_estimator.error_norm.type(var) == H1_Z_SEMINORM || 00883 error_estimator.error_norm.type(var) == H2_SEMINORM) 00884 new_error_per_cell[e] += error_estimator.error_norm.weight_sq(var) * element_error; 00885 else 00886 libmesh_error(); 00887 } // End (re) loop over patch elements 00888 00889 } // end variables loop 00890 00891 // Now that we have the contributions from each variable, 00892 // we have take square roots of the entries we 00893 // added to error_per_cell to get an error norm 00894 // If we are reusing patches, once again reuse the current patch to loop 00895 // over all elements in the current patch, otherwise build a new 00896 // patch containing just the current element and loop over it 00897 Patch::const_iterator patch_re_it; 00898 Patch::const_iterator patch_re_end; 00899 00900 // Build a new patch if necessary 00901 Patch current_elem_patch; 00902 00903 if(this->error_estimator.patch_reuse) 00904 { 00905 // Just get the iterators from the current patch 00906 patch_re_it = patch.begin(); 00907 patch_re_end = patch.end(); 00908 } 00909 else 00910 { 00911 // Use a target patch size of just 0, this will contain 00912 // just the current element. 00913 current_elem_patch.build_around_element (elem, 0, 00914 error_estimator.patch_growth_strategy); 00915 00916 // Get the iterators from this newly constructed patch 00917 patch_re_it = current_elem_patch.begin(); 00918 patch_re_end = current_elem_patch.end(); 00919 } 00920 00921 // Loop over every element in the patch we just constructed 00922 for (unsigned int i = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++i) 00923 { 00924 // The pth element in the patch 00925 const Elem* e_p = *patch_re_it; 00926 00927 // We'll need an index into the error vector 00928 const dof_id_type e_p_id = e_p->id(); 00929 00930 // Update the error_per_cell vector for this element 00931 if (error_estimator.error_norm.type(0) == L2 || 00932 error_estimator.error_norm.type(0) == H1_SEMINORM || 00933 error_estimator.error_norm.type(0) == H1_X_SEMINORM || 00934 error_estimator.error_norm.type(0) == H1_Y_SEMINORM || 00935 error_estimator.error_norm.type(0) == H1_Z_SEMINORM || 00936 error_estimator.error_norm.type(0) == H2_SEMINORM) 00937 { 00938 Threads::spin_mutex::scoped_lock acquire(Threads::spin_mtx); 00939 if (!error_per_cell[e_p_id]) 00940 error_per_cell[e_p_id] = 00941 static_cast<ErrorVectorReal>(std::sqrt(new_error_per_cell[i])); 00942 } 00943 else 00944 { 00945 libmesh_assert (error_estimator.error_norm.type(0) == L_INF || 00946 error_estimator.error_norm.type(0) == W1_INF_SEMINORM || 00947 error_estimator.error_norm.type(0) == W2_INF_SEMINORM); 00948 Threads::spin_mutex::scoped_lock acquire(Threads::spin_mtx); 00949 if (!error_per_cell[e_p_id]) 00950 error_per_cell[e_p_id] = 00951 static_cast<ErrorVectorReal>(new_error_per_cell[i]); 00952 } 00953 00954 } // End loop over every element in patch 00955 00956 } // end element loop 00957 00958 } // End () operator definition
Member Data Documentation
const PatchRecoveryErrorEstimator& libMesh::PatchRecoveryErrorEstimator::EstimateError::error_estimator [private] |
Definition at line 135 of file patch_recovery_error_estimator.h.
Referenced by operator()().
Definition at line 136 of file patch_recovery_error_estimator.h.
Referenced by operator()().
const System& libMesh::PatchRecoveryErrorEstimator::EstimateError::system [private] |
Function to set the boolean patch_reuse in case the user wants to change the default behaviour of patch_recovery_error_estimator
Definition at line 134 of file patch_recovery_error_estimator.h.
Referenced by operator()().
The documentation for this class was generated from the following files:
Site Created By: libMesh Developers
Last modified: February 05 2013 19:55:33 UTC
Hosted By: