21 #ifndef __NITSCHE_PERIODIC_METHOD_HPP__
22 #define __NITSCHE_PERIODIC_METHOD_HPP__
54 int gg,
int side,EntityType
type
62 typedef multi_index_container<
71 member<MultiIndexData,int,&MultiIndexData::gG>,
72 member<MultiIndexData,int,&MultiIndexData::sIde>,
73 member<MultiIndexData,EntityType,&MultiIndexData::tYpe>
97 if(data.getIndices().size()==0) PetscFunctionReturn(0);
99 EntityHandle this_face = getNumeredEntFiniteElementPtr()->getEnt();
100 int nb_gauss_pts_on_this_face = data.getN().size1();
101 for(
int ffgg = 0;ffgg<nb_gauss_pts_on_this_face;ffgg++) {
104 pair<CommonData::Container::iterator,bool> p;
111 int nb_shape_fun = data.getN().size2();
112 shape_fun.resize(nb_shape_fun);
113 cblas_dcopy(nb_shape_fun,&data.getN()(ffgg,0),1,&shape_fun[0],1);
114 p_data.
iNdices = data.getIndices();
115 p_data.
dofOrders.resize(data.getFieldDofs().size(),
false);
116 for(
unsigned int dd = 0;
dd<data.getFieldDofs().size();
dd++) {
117 p_data.
dofOrders[
dd] = data.getFieldDofs()[
dd]->getDofOrder();
119 int nb_dofs = data.getFieldData().size();
120 if(
type == MBVERTEX) {
124 for(
int rr = 0;rr<3;rr++) {
126 nb_dofs/3,&data.getFieldData()[rr],3,&data.getN()(ffgg,0),1
135 for(
int rr = 0;rr<3;rr++) {
137 nb_dofs/3,&data.getFieldData()[rr],3,&data.getN()(ffgg,0),1
142 }
catch (
const std::exception& ex) {
144 ss <<
"throw in method: " << ex.what() << endl;
145 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
149 PetscFunctionReturn(0);
168 EntityHandle face = numeredEntFiniteElementPtr->getEnt();
175 gaussPts(2,ffgg) = 0;
178 }
catch (
const std::exception& ex) {
180 ss <<
"throw in method: " << ex.what() << endl;
181 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
183 PetscFunctionReturn(0);
204 if(data.getIndices().size()==0) PetscFunctionReturn(0);
206 EntityHandle this_tet = getNumeredEntFiniteElementPtr()->getEnt();
207 int nb_gauss_pts = data.getN().size1();
208 for(
int ffgg = 0;ffgg<nb_gauss_pts;ffgg++) {
211 pair<CommonData::Container::iterator,bool> p;
218 diff_shape_fun = data.getDiffN(ffgg);
219 p_data.
iNdices = data.getIndices();
220 if(
type == MBVERTEX) {
225 }
catch (
const std::exception& ex) {
227 ss <<
"throw in method: " << ex.what() << endl;
228 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
230 PetscFunctionReturn(0);
241 VolumeElementForcesAndSourcesCore(m_field),
249 gaussPts.resize(4,0,
false);
250 EntityHandle tet = numeredEntFiniteElementPtr->getEnt();
252 for(
int ff = 0;ff<4;ff++) {
257 const double coords_tet[12] = { 0,0,0, 1,0,0, 0,1,0, 0,0,1 };
259 gaussPts.resize(4,nb_face_gauss_pts);
261 for(
int ffgg = 0;ffgg!=nb_face_gauss_pts;ffgg++,vffgg++) {
266 for(
int dd = 0;
dd<3;
dd++) {
268 N[0]*coords_tet[3*dataH1.facesNodes(ff,0)+
dd]+
269 N[1]*coords_tet[3*dataH1.facesNodes(ff,1)+
dd]+
270 N[2]*coords_tet[3*dataH1.facesNodes(ff,2)+
dd];
272 gaussPts(3,vffgg) = 0;
279 }
catch (
const std::exception& ex) {
281 ss <<
"throw in method: " << ex.what() << endl;
282 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
284 PetscFunctionReturn(0);
309 tRee(&m_field.get_moab()) {
322 PetscFunctionReturn(0);
338 const double tol = 1e-18;
351 for(
int ff = 0;ff<4;ff++) {
356 3,ublas::shallow_array_adaptor<double>(
360 vector<EntityHandle> triangles_out;
361 vector<double> distance_out;
362 rval =
tRee.ray_intersect_triangles(
367 vector<EntityHandle> triangles_out_on_other_side;
368 vector<double> distance_out_on_other_side;
370 vector<EntityHandle>::iterator hit = triangles_out.begin();
371 vector<double>::iterator dit = distance_out.begin();
372 for(;dit!=distance_out.end();dit++,hit++) {
380 *dit *= cblas_ddot(3,normal,1,&ray[0],1)/cblas_dnrm2(3,normal,1);
382 triangles_out_on_other_side.push_back(*hit);
383 distance_out_on_other_side.push_back(*dit);
387 vector<double>::iterator max_dit = max_element(
388 distance_out_on_other_side.begin(),distance_out_on_other_side.end()
391 std::distance(distance_out_on_other_side.begin(),max_dit)
396 local_coords.resize(2);
408 d_local_coords.resize(2,
false);
410 double N[number_nodes];
411 double diffNTRI[number_nodes*2];
412 double nrm2,nrm2_res;
413 const double eps = 1e-14;
414 const int max_it = 10;
417 local_coords.clear();
419 if(number_nodes == 3) {
421 }
else if(number_nodes == 6) {
427 for(
int dd = 0;
dd<3;
dd++) {
430 nrm2_res = cblas_dnrm2(3,&res[0],1);
431 if(number_nodes == 3) {
438 for(
int dd = 0;
dd<3;
dd++) {
443 ATA = prod(trans(
A),
A);
445 ATres = prod(trans(
A),res);
446 double det = (ATA(0,0)*ATA(1,1)-ATA(0,1)*ATA(1,0));
447 d_local_coords[0] = -(ATA(1,1)*ATres[0]-ATA(0,1)*ATres[1])/det;
448 d_local_coords[1] = -(ATA(0,0)*ATres[1]-ATA(1,0)*ATres[0])/det;
449 local_coords += d_local_coords;
450 nrm2 = cblas_dnrm2(2,&d_local_coords[0],1);
452 cerr <<
"Warrning: max_it reached in PeriodicNitscheConstrains:: ... ::doAdditionalJobWhenGuassPtsAreCalulated\n";
453 cerr <<
"nrm2 " << nrm2 <<
" nrm2_res " << nrm2_res << endl;
454 cerr <<
"res " << res << endl;
455 cerr <<
"other_coords " << other_coords << endl;
458 }
while(nrm2>
eps || nrm2_res>
eps);
459 if(local_coords[0]<-
eps) {
462 if(local_coords[1]<-
eps) {
465 if(local_coords[0]>1+
eps) {
468 if(local_coords[1]>1+
eps) {
471 if(local_coords[0]+local_coords[1]>1+
eps) {
475 "local_coords[0]+local_coords[1]>1",
476 local_coords[0]+local_coords[1]
492 }
catch (
const std::exception& ex) {
494 ss <<
"throw in method: " << ex.what() << endl;
495 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
505 NumeredEntFiniteElement_multiIndex::index<Composite_Name_And_Ent_mi_tag>::type::iterator it;
508 .get<Composite_Name_And_Ent_mi_tag>()
513 .get<Composite_Name_And_Ent_mi_tag>()
521 boost::shared_ptr<const NumeredEntFiniteElement> face_ptr = *it;
532 ierr =
mField.get_adjacencies_equality(mit->first,3,tets); CHKERRQ(
ierr);
541 .get<Composite_Name_And_Ent_mi_tag>()
546 .get<Composite_Name_And_Ent_mi_tag>()
553 boost::shared_ptr<const NumeredEntFiniteElement> tet_ptr = *it;
565 }
catch (
const std::exception& ex) {
567 ss <<
"throw in method: " << ex.what() << endl;
568 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
571 PetscFunctionReturn(0);
589 bool field_disp =
true
614 mat_d.resize(3,6,
false);
616 mat_d(0,0) = coords[0]; mat_d(0,3) = coords[1]*0.5; mat_d(0,5) = coords[2]*0.5;
617 mat_d(1,1) = coords[1]; mat_d(1,3) = coords[0]*0.5; mat_d(1,4) = coords[2]*0.5;
618 mat_d(2,2) = coords[2]; mat_d(2,4) = coords[1]*0.5; mat_d(2,5) = coords[0]*0.5;
619 PetscFunctionReturn(0);
628 PetscFunctionReturn(0);
631 if(row_data.getIndices().size()==0) PetscFunctionReturn(0);
632 int nb_dofs = row_data.getIndices().size();
640 for(
int ff = 0;ff<4;ff++) {
644 for(
int ii = 0;ii<6;ii++) {
645 nF0[ii].resize(nb_dofs,
false);
647 nF1[ii].resize(nb_dofs,
false);
650 for(
int fgg = 0;fgg<nb_face_gauss_pts;fgg++,gg++) {
656 3,ublas::shallow_array_adaptor<double>(
660 double area = cblas_dnrm2(3,&normal[0],1);
673 const double err = 1e-10;
674 double delta = inner_prod(other_x,normal)+inner_prod(x,normal);
675 if(fabs(
delta)>err) {
677 cerr << other_x << endl;
683 for(
int ii = 0;ii<6;ii++) {
688 for(
int dd1 = 0;dd1<nb_dofs/3;dd1++) {
689 double n_val = row_data.getN(gg)[dd1];
690 for(
int dd2 = 0;dd2<3;dd2++) {
695 for(
int dd = 0;
dd<nb_dofs;
dd++) {
700 dot -= (1-
eps)*cblas_ddot(
708 for(
int ii = 0;ii<6;ii++) {
713 &row_data.getIndices()[0],
721 &row_data.getIndices()[0],
729 }
catch (
const std::exception& ex) {
731 ss <<
"throw in method: " << ex.what() << endl;
732 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
735 PetscFunctionReturn(0);
751 bool field_disp =
true
775 PetscFunctionReturn(0);
777 if(row_data.getIndices().size()==0) PetscFunctionReturn(0);
778 int nb_dofs_row = row_data.getIndices().size();
786 for(
int ff = 0;ff<4;ff++) {
791 P[ff].resize(nb_face_gauss_pts);
792 for(
int fgg = 0;fgg<nb_face_gauss_pts;fgg++) {
793 P[ff][fgg].resize(3,3);
795 for(
int dd = 0;
dd<3;
dd++) {
796 P[ff][fgg](
dd,
dd) = 1;
800 for(
int fgg = 0;fgg<nb_face_gauss_pts;fgg++,gg++) {
808 double area = cblas_dnrm2(3,&normal[0],1);
809 const double alpha = 0.5;
811 CommonData::Container::nth_index<0>::type::iterator it,hi_it;
814 for(;it!=hi_it;it++) {
816 const ublas::vector<int> &indices = it->iNdices;
819 int nb_dofs_col = indices.size();
820 if(indices.size()==0) {
823 kMatrix0.resize(nb_dofs_row,nb_dofs_col,
false);
825 for(
int dd1 = 0;dd1<nb_dofs_row/3;dd1++) {
826 double n_row = row_data.getN()(gg,dd1);
827 for(
int dd2 = 0;dd2<nb_dofs_col/3;dd2++) {
828 double n_col = shape[dd2];
829 for(
int dd3 = 0;dd3<3;dd3++) {
830 kMatrix0(3*dd1+dd3,3*dd2+dd3) -= (1-
eps)*val*n_row*n_col*area;
834 kMatrix1.resize(nb_dofs_row,nb_dofs_col,
false);
836 for(
int dd1 = 0;dd1<nb_dofs_row;dd1++) {
837 for(
int dd2 = 0;dd2<nb_dofs_col/3;dd2++) {
838 double n_col = shape[dd2];
839 for(
int dd3 = 0;dd3<3;dd3++) {
849 &row_data.getIndices()[0],
859 &row_data.getIndices()[0],
866 kMatrix1.resize(nb_dofs_col,nb_dofs_row);
868 for(
int dd1 = 0;dd1<nb_dofs_col/3;dd1++) {
869 double n_col = shape[dd1];
870 for(
int dd2 = 0;dd2<3;dd2++) {
871 for(
int dd3 = 0;dd3<nb_dofs_row;dd3++) {
883 &row_data.getIndices()[0],
895 for(;it!=hi_it;it++) {
897 int nb = it->iNdices.size();
898 jac.resize(9,nb,
false);
902 for(
int dd = 0;
dd<nb/3;
dd++) {
903 for(
int rr = 0;rr<3;rr++) {
904 for(
int ii = 0;ii<9;ii++) {
905 for(
int jj = 0;jj<3;jj++) {
906 jac(ii,3*
dd+rr) += jac_stress(ii,3*rr+jj)*diff_n(
dd,jj);
911 trac.resize(3,jac.size2());
913 for(
unsigned int dd2 = 0;dd2<jac.size2();dd2++) {
914 for(
int nn = 0;nn<3;nn++) {
915 trac(nn,dd2) = -cblas_ddot(3,&jac(3*nn,dd2),jac.size2(),&normal[0],1);
921 for(
int dd1 = 0;dd1<nb_dofs_row/3;dd1++) {
922 double n_row = row_data.getN()(gg,dd1);
923 for(
int dd2 = 0;dd2<3;dd2++) {
924 for(
int dd3 = 0;dd3<nb;dd3++) {
925 kMatrixOff(3*dd1+dd2,dd3) += (1-
eps)*val*n_row*trac(dd2,dd3);
934 &row_data.getIndices()[0],
946 }
catch (
const std::exception& ex) {
948 ss <<
"throw in method: " << ex.what() << endl;
949 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
952 PetscFunctionReturn(0);
965 Vec stress_homo_ghost_vector
982 if(row_type != MBVERTEX) PetscFunctionReturn(0);
984 PetscFunctionReturn(0);
992 for(
unsigned int gg = 0;gg<row_data.getN().size1();gg++) {
995 if(getHoGaussPtsDetJac().size()>0) {
996 val *= getHoGaussPtsDetJac()[gg];
1001 homoStress[3] += (stress(0,1)+stress(1,0))*0.5*val;
1002 homoStress[4] += (stress(1,2)+stress(2,1))*0.5*val;
1003 homoStress[5] += (stress(0,2)+stress(2,0))*0.5*val;
1006 const int indices_6[6] = {0, 1, 2, 3, 4, 5};
1011 }
catch (
const std::exception& ex) {
1013 ss <<
"throw in method: " << ex.what() << endl;
1014 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1017 PetscFunctionReturn(0);
1032 Vec stress_homo_ghost_vector,
1038 nitsche_common_data,
1061 if(row_type != MBVERTEX) PetscFunctionReturn(0);
1063 PetscFunctionReturn(0);
1072 for(
int ff = 0;ff<4;ff++) {
1075 for(
int fgg = 0;fgg<nb_face_gauss_pts;fgg++,gg++) {
1079 3,ublas::shallow_array_adaptor<double>(
1085 3,ublas::shallow_array_adaptor<double>(
1091 noalias(
tRaction) = prod(normal,stress);
1102 const int indices_6[6] = {0, 1, 2, 3, 4, 5};
1112 }
catch (
const std::exception& ex) {
1114 ss <<
"throw in method: " << ex.what() << endl;
1115 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1118 PetscFunctionReturn(0);
1126 #endif //__NITSCHE_PERIODIC_METHOD_HPP__