33 #ifndef __H__UG__LIB_DISC__SPATIAL_DISC__CONSTRAINTS__LAGRANGE_DIRICHLET_BOUNDARY_IMPL__
34 #define __H__UG__LIB_DISC__SPATIAL_DISC__CONSTRAINTS__LAGRANGE_DIRICHLET_BOUNDARY_IMPL__
51 template <
typename TDomain,
typename TAlgebra>
55 base_type::set_approximation_space(approxSpace);
56 m_spApproxSpace = approxSpace;
57 m_spDomain = approxSpace->domain();
58 m_aaPos = m_spDomain->position_accessor();
61 template <
typename TDomain,
typename TAlgebra>
65 m_vBNDNumberData.clear();
66 m_vNumberData.clear();
67 m_vConstNumberData.clear();
68 m_vVectorData.clear();
71 template <
typename TDomain,
typename TAlgebra>
78 template <
typename TDomain,
typename TAlgebra>
83 for(
size_t i = 0; i < Fcts.size(); ++i){
84 if(i > 0)
function.append(
",");
85 function.append(Fcts[i]);
88 for(
size_t i = 0; i < Subsets.size(); ++i){
89 if(i > 0) subsets.append(
",");
90 subsets.append(Subsets[i]);
93 add(
func,
function.c_str(), subsets.c_str());
96 template <
typename TDomain,
typename TAlgebra>
103 template <
typename TDomain,
typename TAlgebra>
107 std::string
function;
108 for(
size_t i = 0; i < Fcts.size(); ++i){
109 if(i > 0)
function.append(
",");
110 function.append(Fcts[i]);
113 for(
size_t i = 0; i < Subsets.size(); ++i){
114 if(i > 0) subsets.append(
",");
115 subsets.append(Subsets[i]);
118 add(
func,
function.c_str(), subsets.c_str());
121 template <
typename TDomain,
typename TAlgebra>
123 add(
number value,
const char*
function,
const char* subsets)
125 m_vConstNumberData.push_back(
ConstNumberData(value,
function, subsets));
128 template <
typename TDomain,
typename TAlgebra>
130 add(
number value,
const std::vector<std::string>& Fcts,
const std::vector<std::string>& Subsets)
132 std::string
function;
133 for(
size_t i = 0; i < Fcts.size(); ++i){
134 if(i > 0)
function.append(
",");
135 function.append(Fcts[i]);
138 for(
size_t i = 0; i < Subsets.size(); ++i){
139 if(i > 0) subsets.append(
",");
140 subsets.append(Subsets[i]);
143 add(value,
function.c_str(), subsets.c_str());
146 template <
typename TDomain,
typename TAlgebra>
153 template <
typename TDomain,
typename TAlgebra>
157 std::string
function;
158 for(
size_t i = 0; i < Fcts.size(); ++i){
159 if(i > 0)
function.append(
",");
160 function.append(Fcts[i]);
163 for(
size_t i = 0; i < Subsets.size(); ++i){
164 if(i > 0) subsets.append(
",");
165 subsets.append(Subsets[i]);
168 add(
func,
function.c_str(), subsets.c_str());
172 template <
typename TDomain,
typename TAlgebra>
174 add(
const char*
name,
const char*
function,
const char* subsets)
179 add(sp,
function, subsets);
185 add(sp,
function, subsets);
188 if(LuaUserData<MathVector<dim>,
dim>::check_callback_returns(
name)){
191 add(sp,
function, subsets);
197 UG_THROW(
"LagrangeDirichlet::add: Lua-Callback with name '"<<
name<<
198 "' does not exist.");
201 UG_THROW(
"LagrangeDirichlet::add: Cannot find matching callback "
202 "signature. Use one of:\n"
203 "a) Number - Callback\n"
205 "b) Conditional Number - Callback\n"
207 "c) "<<
dim<<
"d Vector - Callback\n"
208 << (LuaUserData<MathVector<dim>,
dim>::signature()));
211 template <
typename TDomain,
typename TAlgebra>
213 add(
const char*
name,
const std::vector<std::string>& Fcts,
const std::vector<std::string>& Subsets)
215 std::string
function;
216 for(
size_t i = 0; i < Fcts.size(); ++i){
217 if(i > 0)
function.append(
",");
218 function.append(Fcts[i]);
221 for(
size_t i = 0; i < Subsets.size(); ++i){
222 if(i > 0) subsets.append(
",");
223 subsets.append(Subsets[i]);
226 add(
name,
function.c_str(), subsets.c_str());
229 template <
typename TDomain,
typename TAlgebra>
231 add(LuaFunctionHandle fct,
const char*
function,
const char* subsets)
235 make_sp(
new LuaUserData<number, dim>(fct));
236 add(sp,
function, subsets);
241 make_sp(
new LuaUserData<number, dim, bool>(fct));
242 add(sp,
function, subsets);
245 if(LuaUserData<MathVector<dim>,
dim>::check_callback_returns(fct)){
247 make_sp(
new LuaUserData<MathVector<dim>,
dim>(fct));
248 add(sp,
function, subsets);
253 UG_THROW(
"LagrangeDirichlet::add: Cannot find matching callback "
254 "signature. Use one of:\n"
255 "a) Number - Callback\n"
257 "b) Conditional Number - Callback\n"
259 "c) "<<
dim<<
"d Vector - Callback\n"
260 << (LuaUserData<MathVector<dim>,
dim>::signature()));
263 template <
typename TDomain,
typename TAlgebra>
265 add(LuaFunctionHandle fct,
const std::vector<std::string>& Fcts,
const std::vector<std::string>& Subsets)
267 std::string
function;
268 for(
size_t i = 0; i < Fcts.size(); ++i){
269 if(i > 0)
function.append(
",");
270 function.append(Fcts[i]);
273 for(
size_t i = 0; i < Subsets.size(); ++i){
274 if(i > 0) subsets.append(
",");
275 subsets.append(Subsets[i]);
278 add(fct,
function.c_str(), subsets.c_str());
282 template <
typename TDomain,
typename TAlgebra>
284 add(
const char* functions,
const char* subsets)
286 m_vOldNumberData.push_back(
OldNumberData(functions, subsets));
289 template <
typename TDomain,
typename TAlgebra>
291 add(
const std::vector<std::string>& Fcts,
const std::vector<std::string>& Subsets)
293 std::string
function;
294 for(
size_t i = 0; i < Fcts.size(); ++i){
295 if(i > 0)
function.append(
",");
296 function.append(Fcts[i]);
299 for(
size_t i = 0; i < Subsets.size(); ++i){
300 if(i > 0) subsets.append(
",");
301 subsets.append(Subsets[i]);
304 add(
function.c_str(), subsets.c_str());
307 template <
typename TDomain,
typename TAlgebra>
312 if(functionGroup.
size() != numFct)
313 UG_THROW(
"DirichletBoundary:extract_data:"
314 " Only "<<numFct<<
" function(s) allowed in specification of a"
315 " Dirichlet Value, but the following functions given:"
322 for(
size_t si = 0; si < subsetGroup.
size(); ++si)
325 const int subsetIndex = subsetGroup[si];
328 if(subsetIndex < 0 || subsetIndex >= pSH->num_subsets())
329 UG_THROW(
"DirichletBoundary:extract_data:"
330 " Invalid Subset Index " << subsetIndex <<
". (Valid is"
331 " 0, .. , " << pSH->num_subsets() <<
").");
334 for(
size_t i=0; i < functionGroup.
size(); ++i)
336 const size_t fct = functionGroup[i];
339 if(fct >= m_spApproxSpace->num_fct())
340 UG_THROW(
"DirichletBoundary:extract_data:"
341 " Function "<< fct <<
" does not exist in pattern.");
344 if(!m_spApproxSpace->is_def_in_subset(fct, subsetIndex))
345 UG_THROW(
"DirichletBoundary:extract_data:"
346 " Function "<<fct<<
" not defined on subset "<<subsetIndex);
353 template <
typename TDomain,
typename TAlgebra>
354 template <
typename TUserData,
typename TScheduledUserData>
356 extract_data(std::map<
int, std::vector<TUserData*> >& mvUserDataBndSegment,
357 std::vector<TScheduledUserData>& vUserData)
360 mvUserDataBndSegment.clear();
362 for(
size_t i = 0; i < vUserData.size(); ++i)
367 if (! m_bInvertSubsetSelection)
368 vUserData[i].ssGrp = m_spApproxSpace->subset_grp_by_name
369 (vUserData[i].ssName.c_str());
371 vUserData[i].ssGrp = m_spApproxSpace->all_subsets_grp_except_for
372 (vUserData[i].ssName.c_str());
375 " all contained in ApproximationSpace.");
379 fctGrp = m_spApproxSpace->fct_grp_by_name(vUserData[i].fctName.c_str());
381 " all contained in ApproximationSpace.");
384 check_functions_and_subsets(fctGrp, vUserData[i].ssGrp, TUserData::numFct);
387 if(fctGrp.
size() != TUserData::numFct)
388 UG_THROW(
"LagrangeDirichletBoundary: wrong number of fct");
390 for(
size_t fct = 0; fct < TUserData::numFct; ++fct)
392 vUserData[i].fct[fct] = fctGrp[fct];
396 for(
size_t si = 0; si < vUserData[i].ssGrp.size(); ++si)
399 const int subsetIndex = vUserData[i].ssGrp[si];
402 mvUserDataBndSegment[subsetIndex].push_back(&vUserData[i]);
408 template <
typename TDomain,
typename TAlgebra>
413 if(!m_spApproxSpace.valid())
414 UG_THROW(
"DirichletBoundary:extract_data: "
415 " Approximation Space not set.");
417 extract_data(m_mNumberBndSegment, m_vNumberData);
418 extract_data(m_mBNDNumberBndSegment, m_vBNDNumberData);
419 extract_data(m_mConstNumberBndSegment, m_vConstNumberData);
420 extract_data(m_mVectorBndSegment, m_vVectorData);
421 extract_data(m_mOldNumberBndSegment, m_vOldNumberData);
428 template <
typename TDomain,
typename TAlgebra>
435 typename std::map<int, std::vector<CondNumberData*> >::const_iterator iter;
436 for(iter = m_mBNDNumberBndSegment.begin(); iter != m_mBNDNumberBndSegment.end(); ++iter)
438 int si = (*iter).first;
439 const std::vector<CondNumberData*>& userData = (*iter).second;
445 std::vector<DoFIndex> multInd;
458 corner = m_aaPos[vertex];
461 for(
size_t i = 0; i < userData.size(); ++i)
464 if(!(*userData[i])(val, corner, time, si))
continue;
467 const size_t fct = userData[i]->fct[0];
470 if(dd->inner_dof_indices(vertex, fct, multInd) != 1)
473 this->m_spAssTuner->set_dirichlet_row(mat, multInd[0]);
484 template <
typename TDomain,
typename TAlgebra>
492 #ifdef LAGRANGE_DIRICHLET_ADJ_TRANSFER_FIX
493 if (!m_bAdjustTransfers)
495 std::cerr <<
"Avoiding adjust_prolongation" << std::endl;
501 adjust_prolongation<CondNumberData>(m_mBNDNumberBndSegment, P, ddFine, ddCoarse, time);
502 adjust_prolongation<NumberData>(m_mNumberBndSegment, P, ddFine, ddCoarse, time);
503 adjust_prolongation<ConstNumberData>(m_mConstNumberBndSegment, P, ddFine, ddCoarse, time);
505 adjust_prolongation<VectorData>(m_mVectorBndSegment, P, ddFine, ddCoarse, time);
507 adjust_prolongation<OldNumberData>(m_mOldNumberBndSegment, P, ddFine, ddCoarse, time);
510 template <
typename TDomain,
typename TAlgebra>
511 template <
typename TUserData>
520 typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
521 for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
524 const int si = (*iter).first;
527 const std::vector<TUserData*>& vUserData = (*iter).second;
532 if(ddFine->max_dofs(
VERTEX)) adjust_prolongation<RegularVertex, TUserData>(vUserData, si, P, ddFine, ddCoarse, time);
533 if(ddFine->max_dofs(
EDGE)) adjust_prolongation<Edge, TUserData>(vUserData, si, P, ddFine, ddCoarse, time);
534 if(ddFine->max_dofs(
FACE)) adjust_prolongation<Face, TUserData>(vUserData, si, P, ddFine, ddCoarse, time);
535 if(ddFine->max_dofs(
VOLUME)) adjust_prolongation<Volume, TUserData>(vUserData, si, P, ddFine, ddCoarse, time);
538 " While calling 'adapt_prolongation' for TUserData, aborting.");
542 template <
typename TDomain,
typename TAlgebra>
543 template <
typename TBaseElem,
typename TUserData>
552 std::vector<DoFIndex> vFineDoF, vCoarseDoF;
558 std::vector<position_type> vPos;
566 for( ; iter != iterEnd; iter++)
570 GridObject* parent = m_spDomain->grid()->get_parent(elem);
571 if(!parent)
continue;
572 if(!ddCoarse->is_contained(parent))
continue;
575 for(
size_t i = 0; i < vUserData.size(); ++i)
577 for(
size_t f = 0; f < TUserData::numFct; ++f)
580 const size_t fct = vUserData[i]->fct[f];
583 const LFEID& lfeID = ddFine->local_finite_element_id(fct);
586 ddFine->inner_dof_indices(elem, fct, vFineDoF);
587 ddCoarse->inner_dof_indices(parent, fct, vCoarseDoF);
590 if(TUserData::isConditional){
591 InnerDoFPosition<TDomain>(vPos, elem, *m_spDomain, lfeID);
592 UG_ASSERT(vFineDoF.size() == vPos.size(),
"Size mismatch");
596 for(
size_t j = 0; j < vFineDoF.size(); ++j)
599 if(TUserData::isConditional){
600 if(!(*vUserData[i])(val, vPos[j], time, si))
continue;
603 SetRow(P, vFineDoF[j], 0.0);
606 if(vFineDoF.size() > 0){
607 for(
size_t k = 0; k < vCoarseDoF.size(); ++k){
608 DoFRef(P, vFineDoF[0], vCoarseDoF[k]) = 1.0;
616 template <
typename TDomain,
typename TAlgebra>
624 #ifdef LAGRANGE_DIRICHLET_ADJ_TRANSFER_FIX
625 if (!m_bAdjustTransfers)
627 std::cerr <<
"Avoiding adjust_restriction" << std::endl;
633 adjust_restriction<CondNumberData>(m_mBNDNumberBndSegment, R, ddCoarse, ddFine, time);
634 adjust_restriction<NumberData>(m_mNumberBndSegment, R, ddCoarse, ddFine, time);
635 adjust_restriction<ConstNumberData>(m_mConstNumberBndSegment, R, ddCoarse, ddFine, time);
637 adjust_restriction<VectorData>(m_mVectorBndSegment, R, ddCoarse, ddFine, time);
639 adjust_restriction<OldNumberData>(m_mOldNumberBndSegment, R, ddCoarse, ddFine, time);
642 template <
typename TDomain,
typename TAlgebra>
643 template <
typename TUserData>
652 typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
653 for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
656 const int si = (*iter).first;
659 const std::vector<TUserData*>& vUserData = (*iter).second;
664 if(ddFine->max_dofs(
VERTEX)) adjust_restriction<RegularVertex, TUserData>(vUserData, si, R, ddCoarse, ddFine, time);
665 if(ddFine->max_dofs(
EDGE)) adjust_restriction<Edge, TUserData>(vUserData, si, R, ddCoarse, ddFine, time);
666 if(ddFine->max_dofs(
FACE)) adjust_restriction<Face, TUserData>(vUserData, si, R, ddCoarse, ddFine, time);
667 if(ddFine->max_dofs(
VOLUME)) adjust_restriction<Volume, TUserData>(vUserData, si, R, ddCoarse, ddFine, time);
670 " While calling 'adjust_restriction' for TUserData, aborting.");
674 template <
typename TDomain,
typename TAlgebra>
675 template <
typename TBaseElem,
typename TUserData>
684 std::vector<DoFIndex> vFineDoF, vCoarseDoF;
690 std::vector<position_type> vPos;
698 for( ; iter != iterEnd; iter++)
702 GridObject* parent = m_spDomain->grid()->get_parent(elem);
703 if(!parent)
continue;
704 if(!ddCoarse->is_contained(parent))
continue;
707 for(
size_t i = 0; i < vUserData.size(); ++i)
709 for(
size_t f = 0; f < TUserData::numFct; ++f)
712 const size_t fct = vUserData[i]->fct[f];
715 const LFEID& lfeID = ddFine->local_finite_element_id(fct);
718 ddFine->inner_dof_indices(elem, fct, vFineDoF);
719 ddCoarse->inner_dof_indices(parent, fct, vCoarseDoF);
722 if(TUserData::isConditional){
723 InnerDoFPosition<TDomain>(vPos, parent, *m_spDomain, lfeID);
724 UG_ASSERT(vCoarseDoF.size() == vPos.size(),
"Size mismatch");
728 for(
size_t j = 0; j < vCoarseDoF.size(); ++j)
731 if(TUserData::isConditional){
732 if(!(*vUserData[i])(val, vPos[j], time, si))
continue;
735 SetRow(R, vCoarseDoF[j], 0.0);
738 if(vFineDoF.size() > 0){
739 for(
size_t k = 0; k < vCoarseDoF.size(); ++k){
740 DoFRef(R, vCoarseDoF[k], vFineDoF[0]) = 1.0;
752 template <
typename TDomain,
typename TAlgebra>
761 adjust_jacobian<CondNumberData>(m_mBNDNumberBndSegment, J, u, dd, time);
762 adjust_jacobian<NumberData>(m_mNumberBndSegment, J, u, dd, time);
763 adjust_jacobian<ConstNumberData>(m_mConstNumberBndSegment, J, u, dd, time);
765 adjust_jacobian<VectorData>(m_mVectorBndSegment, J, u, dd, time);
767 adjust_jacobian<OldNumberData>(m_mOldNumberBndSegment, J, u, dd, time);
771 template <
typename TDomain,
typename TAlgebra>
772 template <
typename TUserData>
774 adjust_jacobian(
const std::map<
int, std::vector<TUserData*> >& mvUserData,
779 typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
780 for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
783 const int si = (*iter).first;
786 const std::vector<TUserData*>& vUserData = (*iter).second;
792 adjust_jacobian<RegularVertex, TUserData>(vUserData, si, J, u, dd, time);
793 if(dd->max_dofs(
EDGE))
794 adjust_jacobian<Edge, TUserData>(vUserData, si, J, u, dd, time);
795 if(dd->max_dofs(
FACE))
796 adjust_jacobian<Face, TUserData>(vUserData, si, J, u, dd, time);
798 adjust_jacobian<Volume, TUserData>(vUserData, si, J, u, dd, time);
801 " While calling 'adapt_jacobian' for TUserData, aborting.");
805 template <
typename TDomain,
typename TAlgebra>
806 template <
typename TBaseElem,
typename TUserData>
813 std::vector<DoFIndex> multInd;
819 std::vector<position_type> vPos;
822 std::set<size_t> dirichletDoFIndices;
831 for( ; iter != iterEnd; iter++)
837 for(
size_t i = 0; i < vUserData.size(); ++i)
839 for(
size_t f = 0; f < TUserData::numFct; ++f)
842 const size_t fct = vUserData[i]->fct[f];
845 const LFEID& lfeID = dd->local_finite_element_id(fct);
848 dd->inner_dof_indices(elem, fct, multInd);
851 if(TUserData::isConditional){
852 InnerDoFPosition<TDomain>(vPos, elem, *m_spDomain, lfeID);
853 UG_ASSERT(multInd.size() == vPos.size(),
"Size mismatch");
857 for(
size_t j = 0; j < multInd.size(); ++j)
860 if(TUserData::isConditional){
861 if(!(*vUserData[i])(val, vPos[j], time, si))
continue;
864 this->m_spAssTuner->set_dirichlet_row(J, multInd[j]);
865 if(m_bDirichletColumns)
866 dirichletDoFIndices.insert(multInd[j][0]);
873 if(m_bDirichletColumns){
877 size_t nr = J.num_rows();
882 typename std::set<size_t>::iterator currentDIndex;
884 for(
size_t i = 0; i<nr; i++)
886 for(
typename matrix_type::row_iterator it = J.begin_row(i); it!=J.end_row(i); ++it){
891 currentDIndex = dirichletDoFIndices.find(it.index());
894 if(currentDIndex != dirichletDoFIndices.end()){
918 template <
typename TDomain,
typename TAlgebra>
923 const std::vector<number>* vScaleMass,
924 const std::vector<number>* vScaleStiff)
928 adjust_defect<CondNumberData>(m_mBNDNumberBndSegment, d, u, dd, time);
929 adjust_defect<NumberData>(m_mNumberBndSegment, d, u, dd, time);
930 adjust_defect<ConstNumberData>(m_mConstNumberBndSegment, d, u, dd, time);
932 adjust_defect<VectorData>(m_mVectorBndSegment, d, u, dd, time);
934 adjust_defect<OldNumberData>(m_mOldNumberBndSegment, d, u, dd, time);
938 template <
typename TDomain,
typename TAlgebra>
939 template <
typename TUserData>
941 adjust_defect(
const std::map<
int, std::vector<TUserData*> >& mvUserData,
946 typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
947 for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
950 const int si = (*iter).first;
953 const std::vector<TUserData*>& vUserData = (*iter).second;
959 adjust_defect<RegularVertex, TUserData>(vUserData, si, d, u, dd, time);
960 if(dd->max_dofs(
EDGE))
961 adjust_defect<Edge, TUserData>(vUserData, si, d, u, dd, time);
962 if(dd->max_dofs(
FACE))
963 adjust_defect<Face, TUserData>(vUserData, si, d, u, dd, time);
965 adjust_defect<Volume, TUserData>(vUserData, si, d, u, dd, time);
968 " While calling 'adjust_defect' for TUserData, aborting.");
972 template <
typename TDomain,
typename TAlgebra>
973 template <
typename TBaseElem,
typename TUserData>
975 adjust_defect(
const std::vector<TUserData*>& vUserData,
int si,
980 std::vector<DoFIndex> multInd;
986 std::vector<position_type> vPos;
994 for( ; iter != iterEnd; iter++)
1000 for(
size_t i = 0; i < vUserData.size(); ++i)
1002 for(
size_t f = 0; f < TUserData::numFct; ++f)
1005 const size_t fct = vUserData[i]->fct[f];
1008 const LFEID& lfeID = dd->local_finite_element_id(fct);
1011 dd->inner_dof_indices(elem, fct, multInd);
1014 if(TUserData::isConditional){
1015 InnerDoFPosition<TDomain>(vPos, elem, *m_spDomain, lfeID);
1016 UG_ASSERT(multInd.size() == vPos.size(),
"Size mismatch. (multInd.size()="<<
1017 multInd.size()<<
", vPos.size()="<<vPos.size()<<
")");
1021 for(
size_t j = 0; j < multInd.size(); ++j)
1024 if(TUserData::isConditional){
1025 if(!(*vUserData[i])(val, vPos[j], time, si))
continue;
1029 this->m_spAssTuner->set_dirichlet_val(d, multInd[j], 0.0);
1040 template <
typename TDomain,
typename TAlgebra>
1046 adjust_solution<CondNumberData>(m_mBNDNumberBndSegment, u, dd, time);
1047 adjust_solution<NumberData>(m_mNumberBndSegment, u, dd, time);
1048 adjust_solution<ConstNumberData>(m_mConstNumberBndSegment, u, dd, time);
1050 adjust_solution<VectorData>(m_mVectorBndSegment, u, dd, time);
1052 adjust_solution<OldNumberData>(m_mOldNumberBndSegment, u, dd, time);
1056 template <
typename TDomain,
typename TAlgebra>
1057 template <
typename TUserData>
1059 adjust_solution(
const std::map<
int, std::vector<TUserData*> >& mvUserData,
1063 typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
1064 for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
1067 const int si = (*iter).first;
1070 const std::vector<TUserData*>& vUserData = (*iter).second;
1076 adjust_solution<RegularVertex, TUserData>(vUserData, si, u, dd, time);
1077 if(dd->max_dofs(
EDGE))
1078 adjust_solution<Edge, TUserData>(vUserData, si, u, dd, time);
1079 if(dd->max_dofs(
FACE))
1080 adjust_solution<Face, TUserData>(vUserData, si, u, dd, time);
1082 adjust_solution<Volume, TUserData>(vUserData, si, u, dd, time);
1085 " While calling 'adjust_solution' for TUserData, aborting.");
1089 template <
typename TDomain,
typename TAlgebra>
1090 template <
typename TBaseElem,
typename TUserData>
1096 if (! TUserData::setSolValue)
1100 std::vector<DoFIndex> multInd;
1106 std::vector<position_type> vPos;
1114 for( ; iter != iterEnd; iter++)
1120 for(
size_t i = 0; i < vUserData.size(); ++i)
1122 for(
size_t f = 0; f < TUserData::numFct; ++f)
1125 const size_t fct = vUserData[i]->fct[f];
1128 const LFEID& lfeID = dd->local_finite_element_id(fct);
1131 InnerDoFPosition<TDomain>(vPos, elem, *m_spDomain, lfeID);
1134 dd->inner_dof_indices(elem, fct, multInd);
1136 UG_ASSERT(multInd.size() == vPos.size(),
"Size mismatch");
1139 for(
size_t j = 0; j < multInd.size(); ++j)
1142 if(!(*vUserData[i])(val, vPos[j], time, si))
continue;
1144 this->m_spAssTuner->set_dirichlet_val(u, multInd[j], val[f]);
1157 template <
typename TDomain,
typename TAlgebra>
1168 adjust_correction<CondNumberData>(m_mBNDNumberBndSegment, c, dd, time);
1169 adjust_correction<NumberData>(m_mNumberBndSegment, c, dd, time);
1170 adjust_correction<ConstNumberData>(m_mConstNumberBndSegment, c, dd, time);
1172 adjust_correction<VectorData>(m_mVectorBndSegment, c, dd, time);
1174 adjust_correction<OldNumberData>(m_mOldNumberBndSegment, c, dd, time);
1177 template <
typename TDomain,
typename TAlgebra>
1178 template <
typename TUserData>
1184 typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
1185 for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
1188 const int si = (*iter).first;
1191 const std::vector<TUserData*>& vUserData = (*iter).second;
1197 adjust_correction<RegularVertex, TUserData>(vUserData, si, c, dd, time);
1198 if(dd->max_dofs(
EDGE))
1199 adjust_correction<Edge, TUserData>(vUserData, si,c, dd, time);
1200 if(dd->max_dofs(
FACE))
1201 adjust_correction<Face, TUserData>(vUserData, si, c, dd, time);
1203 adjust_correction<Volume, TUserData>(vUserData, si, c, dd, time);
1206 " While calling 'adjust_correction' for TUserData, aborting.");
1210 template <
typename TDomain,
typename TAlgebra>
1211 template <
typename TBaseElem,
typename TUserData>
1217 std::vector<DoFIndex> multInd;
1223 std::vector<position_type> vPos;
1231 for( ; iter != iterEnd; iter++)
1237 for(
size_t i = 0; i < vUserData.size(); ++i)
1239 for(
size_t f = 0; f < TUserData::numFct; ++f)
1242 const size_t fct = vUserData[i]->fct[f];
1245 const LFEID& lfeID = dd->local_finite_element_id(fct);
1248 InnerDoFPosition<TDomain>(vPos, elem, *m_spDomain, lfeID);
1251 dd->inner_dof_indices(elem, fct, multInd);
1253 UG_ASSERT(multInd.size() == vPos.size(),
"Size mismatch");
1256 for(
size_t j = 0; j < multInd.size(); ++j)
1259 if(!(*vUserData[i])(val, vPos[j], time, si))
continue;
1261 this->m_spAssTuner->set_dirichlet_val(c, multInd[j], 0.0);
1273 template <
typename TDomain,
typename TAlgebra>
1280 adjust_linear<CondNumberData>(m_mBNDNumberBndSegment, A, b, dd, time);
1281 adjust_linear<NumberData>(m_mNumberBndSegment, A, b, dd, time);
1282 adjust_linear<ConstNumberData>(m_mConstNumberBndSegment, A, b, dd, time);
1284 adjust_linear<VectorData>(m_mVectorBndSegment, A, b, dd, time);
1286 adjust_linear<OldNumberData>(m_mOldNumberBndSegment, A, b, dd, time);
1290 template <
typename TDomain,
typename TAlgebra>
1291 template <
typename TUserData>
1293 adjust_linear(
const std::map<
int, std::vector<TUserData*> >& mvUserData,
1298 typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
1299 for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
1302 const int si = (*iter).first;
1305 const std::vector<TUserData*>& vUserData = (*iter).second;
1311 adjust_linear<RegularVertex, TUserData>(vUserData, si, A, b, dd, time);
1312 if(dd->max_dofs(
EDGE))
1313 adjust_linear<Edge, TUserData>(vUserData, si, A, b, dd, time);
1314 if(dd->max_dofs(
FACE))
1315 adjust_linear<Face, TUserData>(vUserData, si, A, b, dd, time);
1317 adjust_linear<Volume, TUserData>(vUserData, si, A, b, dd, time);
1320 " While calling 'adjust_linear' for TUserData, aborting.");
1324 template <
typename TDomain,
typename TAlgebra>
1325 template <
typename TBaseElem,
typename TUserData>
1327 adjust_linear(
const std::vector<TUserData*>& vUserData,
int si,
1332 std::vector<DoFIndex> multInd;
1338 std::set<size_t> dirichletDoFIndices;
1341 std::vector<position_type> vPos;
1349 for( ; iter != iterEnd; iter++)
1355 for(
size_t i = 0; i < vUserData.size(); ++i)
1357 for(
size_t f = 0; f < TUserData::numFct; ++f)
1360 const size_t fct = vUserData[i]->fct[f];
1363 const LFEID& lfeID = dd->local_finite_element_id(fct);
1366 InnerDoFPosition<TDomain>(vPos, elem, *m_spDomain, lfeID);
1369 dd->inner_dof_indices(elem, fct, multInd);
1371 UG_ASSERT(multInd.size() == vPos.size(),
1372 "Mismatch: numInd="<<multInd.size()<<
", numPos="
1373 <<vPos.size()<<
" on "<<elem->reference_object_id());
1376 for(
size_t j = 0; j < multInd.size(); ++j)
1379 if(!(*vUserData[i])(val, vPos[j], time, si))
continue;
1381 this->m_spAssTuner->set_dirichlet_row(A, multInd[j]);
1383 if(m_bDirichletColumns)
1388 "adjust_linear() is not implemented for block matrices and the symmetric case!");
1389 dirichletDoFIndices.insert(multInd[j][0]);
1392 if (TUserData::setSolValue)
1393 this->m_spAssTuner->set_dirichlet_val(b, multInd[j], val[f]);
1401 if(m_bDirichletColumns){
1405 size_t nr = A.num_rows();
1407 typename std::set<size_t>::iterator currentDIndex;
1412 for(
size_t i = 0; i<nr; i++)
1415 if (dirichletDoFIndices.find(i) != dirichletDoFIndices.end())
1418 for(
typename matrix_type::row_iterator it = A.begin_row(i); it!=A.end_row(i); ++it)
1420 currentDIndex = dirichletDoFIndices.find(it.index());
1423 if(currentDIndex != dirichletDoFIndices.end()){
1427 m_dirichletMap[si][i][it.index()] = it.value();
1434 typename std::map<int, std::map<int, value_type> >::iterator itdirichletMap;
1435 typename std::map<int, value_type>::iterator itdirichletRowMap;
1436 for(
size_t i = 0; i<nr; i++)
1439 if ((itdirichletMap = m_dirichletMap[si].
find(i)) == m_dirichletMap[si].end())
1442 for(
typename matrix_type::row_iterator it = m_A->begin_row(i); it!=m_A->end_row(i); ++it){
1445 if ((itdirichletRowMap = itdirichletMap->second.find(it.index())) != itdirichletMap->second.end())
1446 b[i] -= itdirichletRowMap->second*b[it.index()];
1456 template <
typename TDomain,
typename TAlgebra>
1463 adjust_rhs<CondNumberData>(m_mBNDNumberBndSegment, b, u, dd, time);
1464 adjust_rhs<NumberData>(m_mNumberBndSegment, b, u, dd, time);
1465 adjust_rhs<ConstNumberData>(m_mConstNumberBndSegment, b, u, dd, time);
1467 adjust_rhs<VectorData>(m_mVectorBndSegment, b, u, dd, time);
1469 adjust_rhs<OldNumberData>(m_mOldNumberBndSegment, b, u, dd, time);
1473 template <
typename TDomain,
typename TAlgebra>
1474 template <
typename TUserData>
1476 adjust_rhs(
const std::map<
int, std::vector<TUserData*> >& mvUserData,
1481 typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
1482 for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
1485 const int si = (*iter).first;
1488 const std::vector<TUserData*>& vUserData = (*iter).second;
1494 adjust_rhs<RegularVertex, TUserData>(vUserData, si, b, u, dd, time);
1495 if(dd->max_dofs(
EDGE))
1496 adjust_rhs<Edge, TUserData>(vUserData, si, b, u, dd, time);
1497 if(dd->max_dofs(
FACE))
1498 adjust_rhs<Face, TUserData>(vUserData, si, b, u, dd, time);
1500 adjust_rhs<Volume, TUserData>(vUserData, si, b, u, dd, time);
1503 " While calling 'adjust_rhs' for TUserData, aborting.");
1507 template <
typename TDomain,
typename TAlgebra>
1508 template <
typename TBaseElem,
typename TUserData>
1510 adjust_rhs(
const std::vector<TUserData*>& vUserData,
int si,
1515 std::vector<DoFIndex> multInd;
1521 std::vector<position_type> vPos;
1529 for( ; iter != iterEnd; iter++)
1535 for(
size_t i = 0; i < vUserData.size(); ++i)
1537 for(
size_t f = 0; f < TUserData::numFct; ++f)
1540 const size_t fct = vUserData[i]->fct[f];
1543 const LFEID& lfeID = dd->local_finite_element_id(fct);
1546 InnerDoFPosition<TDomain>(vPos, elem, *m_spDomain, lfeID);
1549 dd->inner_dof_indices(elem, fct, multInd);
1551 UG_ASSERT(multInd.size() == vPos.size(),
"Size mismatch");
1554 for(
size_t j = 0; j < multInd.size(); ++j)
1557 if(!(*vUserData[i])(val, vPos[j], time, si))
continue;
1559 if (TUserData::setSolValue)
1560 this->m_spAssTuner->set_dirichlet_val(b, multInd[j], val[f]);
1562 this->m_spAssTuner->set_dirichlet_val(b, multInd[j],
DoFRef(u, multInd[j]));
1571 if(m_bDirichletColumns){
1572 typename std::map<int, std::map<int, value_type> >::iterator itdirichletMap;
1573 typename std::map<int, value_type>::iterator itdirichletRowMap;
1574 size_t nr = m_A->num_rows();
1575 for(
size_t i = 0; i<nr; i++)
1578 if ((itdirichletMap = m_dirichletMap[si].
find(i)) == m_dirichletMap[si].end())
1581 for(
typename matrix_type::row_iterator it = m_A->begin_row(i); it!=m_A->end_row(i); ++it){
1584 if ((itdirichletRowMap = itdirichletMap->second.find(it.index())) != itdirichletMap->second.end())
1585 b[i] -= itdirichletRowMap->second*b[it.index()];
1596 template <
typename TDomain,
typename TAlgebra>
1605 const std::vector<number>* vScaleMass,
1606 const std::vector<number>* vScaleStiff
1610 if (this->m_spErrEstData.get() == NULL)
1612 UG_THROW(
"No ErrEstData object has been given to this constraint!");
1619 UG_THROW(
"Dynamic cast to MultipleSideAndElemErrEstData failed."
1620 << std::endl <<
"Make sure you handed the correct type of ErrEstData to this discretization.");
1626 adjust_error<CondNumberData>(m_mBNDNumberBndSegment, u, dd, time);
1627 adjust_error<NumberData>(m_mNumberBndSegment, u, dd, time);
1628 adjust_error<ConstNumberData>(m_mConstNumberBndSegment, u, dd, time);
1630 adjust_error<VectorData>(m_mVectorBndSegment, u, dd, time);
1632 adjust_error<OldNumberData>(m_mOldNumberBndSegment, u, dd, time);
1635 template <
typename TDomain,
typename TAlgebra>
1636 template <
typename TUserData>
1640 const std::map<
int, std::vector<TUserData*> >& mvUserData,
1652 typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
1653 for (iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
1656 const int si = (*iter).first;
1659 const std::vector<TUserData*>& vUserData = (*iter).second;
1664 std::vector<DoFIndex> multInd;
1670 std::vector<position_type> vPos;
1674 iter = dd->begin<side_type>(si);
1675 iterEnd = dd->end<side_type>(si);
1679 for( ; iter != iterEnd; iter++)
1682 side_type* elem = *iter;
1688 std::vector<typename TDomain::position_type> vCoCo;
1692 for(
size_t i = 0; i < vUserData.size(); ++i)
1694 for(
size_t f = 0; f < TUserData::numFct; ++f)
1697 const size_t fct = vUserData[i]->fct[f];
1700 LFEID lfeID = dd->local_finite_element_id(fct);
1709 numSideIPs = err_est_data->
get(fct)->num_side_ips(roid);
1710 sideLocIPs = err_est_data->
get(fct)->template side_local_ips<side_type::dim>(roid);
1711 sideGlobIPs = err_est_data->
get(fct)->side_global_ips(elem, &vCoCo[0]);
1713 UG_CATCH_THROW(
"Global integration points for error estimator cannot be determined.");
1716 for (
size_t ip = 0; ip < numSideIPs; ++ip)
1719 if (!(*vUserData[i])(val, sideGlobIPs[ip], time, si))
continue;
1722 if (! TUserData::setSolValue)
1724 (*err_est_data->
get(fct))(elem,ip) = 0;
1731 LocalFiniteElementProvider::get<side_type::dim>(roid, new_lfeID);
1732 std::vector<number> vShape;
1733 rTrialSpace.
shapes(vShape, sideLocIPs[ip]);
1736 std::vector<DoFIndex> ind;
1737 dd->dof_indices(elem, fct, ind);
1741 for (
size_t sh = 0; sh < vShape.size(); ++sh)
1742 uAtIP +=
DoFRef(u, ind[sh]) * vShape[sh];
1745 (*err_est_data->
get(fct))(elem,ip) = val[f] - uAtIP;
1752 " While calling 'adjust_error' for TUserData, aborting.");
location name
Definition: checkpoint_util.lua:128
Definition: smart_pointer.h:296
Definition: smart_pointer.h:108
base class for approximation spaces without type of algebra or dof distribution
Definition: approximation_space.h:279
virtual void adjust_error(const vector_type &u, ConstSmartPtr< DoFDistribution > dd, int type, number time=0.0, ConstSmartPtr< VectorTimeSeries< vector_type > > vSol=SPNULL, const std::vector< number > *vScaleMass=NULL, const std::vector< number > *vScaleStiff=NULL)
adjusts the assembled error estimator values in the attachments according to the constraint
Definition: lagrange_dirichlet_boundary_impl.h:1599
void add(SmartPtr< UserData< number, dim, bool > > func, const char *function, const char *subsets)
adds a lua callback (cond and non-cond)
Definition: lagrange_dirichlet_boundary_impl.h:73
domain_type::position_type position_type
Type of position coordinates (e.g. position_type)
Definition: lagrange_dirichlet_boundary.h:68
void assemble_dirichlet_rows(matrix_type &mat, ConstSmartPtr< DoFDistribution > dd, number time=0.0)
Sets dirichlet rows for all registered dirichlet values.
Definition: lagrange_dirichlet_boundary_impl.h:430
void adjust_linear(matrix_type &A, vector_type &b, ConstSmartPtr< DoFDistribution > dd, int type, number time=0.0)
sets unity rows in A and dirichlet values in right-hand side b
Definition: lagrange_dirichlet_boundary_impl.h:1275
void adjust_defect(vector_type &d, const vector_type &u, ConstSmartPtr< DoFDistribution > dd, int type, number time=0.0, ConstSmartPtr< VectorTimeSeries< vector_type > > vSol=NULL, const std::vector< number > *vScaleMass=NULL, const std::vector< number > *vScaleStiff=NULL)
sets a zero value in the defect for all dirichlet indices
Definition: lagrange_dirichlet_boundary_impl.h:920
void set_approximation_space(SmartPtr< ApproximationSpace< TDomain > > approxSpace)
sets the approximation space to work on
Definition: lagrange_dirichlet_boundary_impl.h:53
virtual void adjust_restriction(matrix_type &R, ConstSmartPtr< DoFDistribution > ddCoarse, ConstSmartPtr< DoFDistribution > ddFine, int type, number time=0.0)
sets constraints in restriction
Definition: lagrange_dirichlet_boundary_impl.h:618
void check_functions_and_subsets(FunctionGroup &functionGroup, SubsetGroup &subsetGroup, size_t numFct) const
Definition: lagrange_dirichlet_boundary_impl.h:309
void adjust_rhs(vector_type &b, const vector_type &u, ConstSmartPtr< DoFDistribution > dd, int type, number time=0.0)
sets the dirichlet value in the right-hand side
Definition: lagrange_dirichlet_boundary_impl.h:1458
void adjust_solution(vector_type &u, ConstSmartPtr< DoFDistribution > dd, int type, number time=0.0)
sets the dirichlet value in the solution for all dirichlet indices
Definition: lagrange_dirichlet_boundary_impl.h:1042
void adjust_jacobian(matrix_type &J, const vector_type &u, ConstSmartPtr< DoFDistribution > dd, int type, number time=0.0, ConstSmartPtr< VectorTimeSeries< vector_type > > vSol=NULL, const number s_a0=1.0)
sets a unity row for all dirichlet indices
Definition: lagrange_dirichlet_boundary_impl.h:754
virtual void adjust_prolongation(matrix_type &P, ConstSmartPtr< DoFDistribution > ddFine, ConstSmartPtr< DoFDistribution > ddCoarse, int type, number time=0.0)
sets constraints in prolongation
Definition: lagrange_dirichlet_boundary_impl.h:486
void clear()
removes all scheduled dirichlet data.
Definition: lagrange_dirichlet_boundary_impl.h:63
void extract_data()
Definition: lagrange_dirichlet_boundary_impl.h:410
virtual void adjust_correction(vector_type &c, ConstSmartPtr< DoFDistribution > dd, int type, number time=0.0)
sets zero to correction
Definition: lagrange_dirichlet_boundary_impl.h:1159
Definition: function_group.h:52
size_t size() const
number of functions in this group
Definition: function_group.h:116
The base class for all geometric objects, such as vertices, edges, faces, volumes,...
Definition: grid_base_objects.h:157
algebra_type::vector_type vector_type
Type of algebra vector.
Definition: constraint_interface.h:72
algebra_type::matrix_type matrix_type
Type of algebra matrix.
Definition: constraint_interface.h:69
Identifier for Local Finite Elements.
Definition: local_finite_element_id.h:98
SpaceType type() const
returns the type of the local finite element
Definition: local_finite_element_id.h:132
int order() const
returns the order of the local finite element
Definition: local_finite_element_id.h:126
int dim() const
returns the space dimension of the local finite element
Definition: local_finite_element_id.h:129
virtual base class for local shape function sets
Definition: local_shape_function_set.h:70
virtual void shapes(shape_type *vShape, const MathVector< dim > &x) const =0
returns all shape functions evaluated at a point
static SmartPtr< LuaUserData< TData, dim, TRet > > create(const std::string &name)
Definition: lua_user_data.h:223
provides data specified in the lua script
Definition: lua_user_data.h:96
static bool check_callback_returns(const char *callName, const bool bThrow=false)
returns true if callback has correct return values
Definition: lua_user_data_impl.h:249
static std::string signature()
returns string of required callback signature
Definition: lua_user_data_impl.h:67
TErrEstData * get(size_t uid)
accessing the underlying error estimator data objects via function id
Definition: err_est_data.h:466
Definition: err_est_data.h:506
SideAndElemErrEstData< TDomain >::side_type side_type
type of the sides (face, edge) and the elems (volume, face)
Definition: err_est_data.h:512
Group of subsets.
Definition: subset_group.h:51
size_t size() const
number of subsets in this group
Definition: subset_group.h:122
Const iterator to traverse the surface of a multi-grid hierarchy.
Definition: surface_view.h:237
time series of solutions and corresponding time point
Definition: solution_time_series.h:59
Base-class for all vertex-types.
Definition: grid_base_objects.h:231
void SetRow(TSparseMatrix &A, size_t i, size_t alpha, number val=0.0)
Definition: sparsematrix_util.h:710
void CollectCornerCoordinates(std::vector< typename TAAPos::ValueType > &vCornerCoordsOut, const TElem &elem, const TAAPos &aaPos, bool clearContainer=true)
returns the corner coordinates of a geometric object
Definition: domain_util_impl.h:75
base_type::TBaseElem TBaseElem
#define UG_ASSERT(expr, msg)
Definition: assert.h:70
#define UG_CATCH_THROW(msg)
Definition: error.h:64
#define UG_THROW(msg)
Definition: error.h:57
#define UG_COND_THROW(cond, msg)
UG_COND_THROW(cond, msg) : performs a UG_THROW(msg) if cond == true.
Definition: error.h:61
double number
Definition: types.h:124
function util LuaCallbackHelper create(func)
IndexLayout::Interface::iterator find(IndexLayout::Interface &interface, size_t i)
Definition: parallel_index_layout.h:77
ReferenceObjectID
these ids are used to identify the shape of a geometric object.
Definition: grid_base_objects.h:74
number & DoFRef(TMatrix &mat, const DoFIndex &iInd, const DoFIndex &jInd)
Definition: multi_index.h:276
@ VOLUME
Definition: grid_base_objects.h:63
@ VERTEX
Definition: grid_base_objects.h:60
@ EDGE
Definition: grid_base_objects.h:61
@ FACE
Definition: grid_base_objects.h:62
bool CheckLuaCallbackName(const char *name)
returns true if callback exists
Definition: lua_user_data.cpp:50
function func(x, y, z, t, si)
SmartPtr< T, FreePolicy > make_sp(T *inst)
returns a SmartPtr for the passed raw pointer
Definition: smart_pointer.h:836
T value_type
Definition: sparsematrix_interface.h:2
grouping for subset and conditional data
Definition: lagrange_dirichlet_boundary.h:363
grouping for subset and constant data
Definition: lagrange_dirichlet_boundary.h:387
grouping for subset and non-conditional data
Definition: lagrange_dirichlet_boundary.h:338
grouping for subset and the data already stored in the solution
Definition: lagrange_dirichlet_boundary.h:435
grouping for subset and non-conditional vector data
Definition: lagrange_dirichlet_boundary.h:411