6 OutPut(
"Example: Hemker_DEVSS" << endl) ;
13 #include <BoundEdge.h>
14 #include <BoundComp.h>
18 #include <MainUtilities.h>
21 #include <MacroCell.h>
22 #include <BoundEdge.h>
23 #include <IsoBoundEdge.h>
25 #include <IsoInterfaceJoint.h>
30 #include <QuadAffin.h>
31 #include <QuadBilinear.h>
32 #include <QuadIsoparametric.h>
33 #include <TriaAffin.h>
34 #include <TriaIsoparametric.h>
45 void ExactU1(
double x,
double y,
double *values)
53 void ExactU2(
double x,
double y,
double *values)
61 void ExactP(
double x,
double y,
double *values)
72 void BoundCondition(
int i,
double t, BoundCond &cond)
88 void U1BoundValue(
int BdComp,
double Param,
double &value)
98 case 3: value=1.5*(1.0 - (0.25*pow(2-4*Param,2))) ;
102 default: cout <<
"wrong boundary part number: " << BdComp << endl;
106 void U2BoundValue(
int BdComp,
double Param,
double &value)
109 if(BdComp>4) cout <<
"wrong boundary part number: " << BdComp << endl;
115 void LinCoeffs(
int n_points,
double *x,
double *y,
116 double **parameters,
double **coeffs)
120 double *coeff, nondim;
127 for(i=0;i<n_points;i++)
142 double &cd,
double &cl)
145 int N_Points,N_Edges,comp,ed_nr;
146 double *weights, *xi, *eta;
147 double X[MaxN_QuadPoints_2D];
148 double Y[MaxN_QuadPoints_2D];
149 double AbsDetjk[MaxN_QuadPoints_2D];
150 int N_LocalUsedElements;
151 FE2D LocalUsedElements[2], CurrentElement;
153 double **OrigFEValues, *Orig;
154 bool SecondDer[2] = { FALSE, FALSE };
157 int *UGlobalNumbers, *UBeginIndex;
158 int *PGlobalNumbers, *PBeginIndex;
159 int *N_BaseFunct, N_Cells;
160 BaseFunct2D BaseFunct, *BaseFuncts;
163 double value, value1, value2, value3;
164 double FEFunctValues[MaxN_BaseFunctions2D];
165 double FEFunctValues1[MaxN_BaseFunctions2D];
166 double FEFunctValues2[MaxN_BaseFunctions2D];
167 double FEFunctValues3[MaxN_BaseFunctions2D];
168 int N_DerivativesU = 3;
169 double *Derivatives[MaxN_BaseFunctions2D];
170 MultiIndex2D NeededDerivatives[3] = { D00, D10, D01 };
180 int N_DOF_Circ, *DOF_Circ;
181 char VString[] =
"v";
199 aux =
new double [MaxN_QuadPoints_2D*10];
200 for(j=0;j<MaxN_QuadPoints_2D;j++)
201 Derivatives[j] = aux + j*10;
205 memset(v,0,N_*SizeOfDouble);
214 for(i=0;i<N_Cells;i++)
218 for(j=0;j<N_Edges;j++)
221 if ((joint->
GetType() == BoundaryEdge)||
222 (joint->
GetType() == IsoBoundEdge))
227 comp=BoundComp->
GetID();
230 FEEle = USpace->
GetFE2D(i,cell);
235 DOF = UGlobalNumbers + UBeginIndex[i];
236 for (k=0;k<N_DOF_Circ;k++)
237 v[DOF[DOF_Circ[k]]] = 1;
251 for(i=0;i<N_Cells;i++)
258 N_LocalUsedElements = 2;
259 LocalUsedElements[0] = USpace->
GetFE2D(i, cell);
260 LocalUsedElements[1] = PSpace->
GetFE2D(i, cell);
265 TFEDatabase2D::GetOrig(N_LocalUsedElements, LocalUsedElements, Coll, cell, SecondDer, N_Points, xi, eta, weights, X, Y, AbsDetjk);
268 CurrentElement = LocalUsedElements[1];
269 BaseFunct = BaseFuncts[CurrentElement];
270 N_ = N_BaseFunct[CurrentElement];
272 DOF = PGlobalNumbers + PBeginIndex[i];
274 FEFunctValues[l] = p[DOF[l]];
278 for(j=0;j<N_Points;j++)
280 Orig = OrigFEValues[j];
283 value += FEFunctValues[l] * Orig[l];
285 Derivatives[j][0] = value;
289 CurrentElement = LocalUsedElements[0];
290 BaseFunct = BaseFuncts[CurrentElement];
291 N_ = N_BaseFunct[CurrentElement];
293 DOF = UGlobalNumbers + UBeginIndex[i];
296 FEFunctValues1[l] = u1[DOF[l]];
297 FEFunctValues2[l] = u2[DOF[l]];
298 FEFunctValues3[l] = v[DOF[l]];
301 for(k=0;k<N_DerivativesU;k++)
304 NeededDerivatives[k]);
305 for(j=0;j<N_Points;j++)
307 Orig = OrigFEValues[j];
313 value1 += FEFunctValues1[l] * Orig[l];
314 value2 += FEFunctValues2[l] * Orig[l];
315 value3 += FEFunctValues3[l] * Orig[l];
317 Derivatives[j][k+1] = value1;
318 Derivatives[j][k+4] = value2;
319 Derivatives[j][k+7] = value3;
324 for(j=0;j<N_Points;j++)
326 Der = Derivatives[j];
329 value1 = nu*(Der[2]*Der[8]+Der[3]*Der[9]);
333 value1 -= Der[0]*Der[8];
335 value2 = (Der[5]*Der[8]+Der[6]*Der[9]);
337 value2 -= Der[0]*Der[9];
339 cd += AbsDetjk[j]*weights[j] * value1;
340 cl += AbsDetjk[j]*weights[j] * value2;
348 delete Derivatives[0];
358 void TriaReMeshGen(
TDomain *&Domain)
360 int j, ID, k, N_G, *PartMarker, *PointNeighb, maxEpV=0;
361 int a, b, len1, len2, Neighb_tmp, BDpart;
362 int i, temp, N_Cells, N_RootCells, CurrVertex, N_Joints, N_Vertices;
363 int N_Interface_Vert, N_Verti, N_Hori, N_SlipBound_Vert, N_BDVertices;
364 int CurrComp, In_Index, *Triangles, Neib[2], CurrNeib;
366 double deviation, hi, x0, y0, x, y, phi1, phi2;
367 double T_a, T_b, C_x, C_y, s, theta;
369 double *Coordinates, *Hole_List;
371 double Xi[4] = {-20., 20., 20.,-20.};
372 double Yi[4] = {-2.,-2., 2., 2.};
378 TVertex **VertexDel, **NewVertices;
386 std::ostringstream opts;
389 Out.pointlist = NULL;
390 Out.pointattributelist = NULL;
391 Out.pointmarkerlist = NULL;
392 Out.trianglelist = NULL;
393 Out.triangleattributelist = NULL;
394 Out.trianglearealist = NULL;
395 Out.neighborlist = NULL;
396 Out.segmentlist = NULL;
397 Out.segmentmarkerlist = NULL;
399 Out.regionlist = NULL;
401 Out.edgemarkerlist = NULL;
404 opts.seekp(std::ios::beg);
434 N_SlipBound_Vert = 2*N_Hori + 2*N_Verti;
436 N_BDVertices = N_Interface_Vert+N_SlipBound_Vert;
437 In.numberofpoints = N_BDVertices;
438 In.pointlist =
new double[2*In.numberofpoints];
439 In.pointmarkerlist =
new int[In.numberofpoints];
440 In.numberofpointattributes = 0;
442 In.numberofsegments = In.numberofpoints;
443 In.segmentlist =
new int[2*In.numberofsegments];
444 In.segmentmarkerlist =
new int[In.numberofsegments];
445 In.numberofregions = 0;
446 In.regionlist = NULL;
448 In.numberofholes = 1;
451 Hole_List =
new double[2* In.numberofholes];
454 In.holelist = Hole_List;
459 hi = (Xi[1] - Xi[0])/(
double)N_Hori;
465 for(i=0;i<N_Hori;i++)
467 x = x0 + (double)i*hi;
468 In.pointlist[2*In_Index] = x;
469 In.pointlist[2*In_Index+1] = y0;
471 In.pointmarkerlist[In_Index] = CurrComp;
472 In.segmentlist[2*In_Index] = In_Index;
473 In.segmentlist[2*In_Index+1] = In_Index+1;
474 In.segmentmarkerlist[In_Index] = CurrComp;
480 hi = (Yi[2] - Yi[1])/(
double)N_Verti;
485 for(i=0;i<N_Verti;i++)
487 y = y0 + (double)i*hi;
488 In.pointlist[2*In_Index] = x0;
489 In.pointlist[2*In_Index+1] = y;
491 In.pointmarkerlist[In_Index] = CurrComp;
492 In.segmentlist[2*In_Index] = In_Index;
493 In.segmentlist[2*In_Index+1] = In_Index+1;
494 In.segmentmarkerlist[In_Index] = CurrComp;
502 hi = (Xi[3] - Xi[2])/(
double)N_Hori;
507 for(i=0;i<N_Hori;i++)
509 x = x0 + (double)i*hi;
510 In.pointlist[2*In_Index] = x;
511 In.pointlist[2*In_Index+1] = y0;
513 In.pointmarkerlist[In_Index] = CurrComp;
514 In.segmentlist[2*In_Index] = In_Index;
515 In.segmentlist[2*In_Index+1] = In_Index+1;
516 In.segmentmarkerlist[In_Index] = CurrComp;
524 hi = (Yi[0] - Yi[3])/(
double)N_Verti;
529 for(i=0;i<N_Verti;i++)
531 y = y0 + (double)i*hi;
532 In.pointlist[2*In_Index] = x0;
533 In.pointlist[2*In_Index+1] = y;
535 In.pointmarkerlist[In_Index] = CurrComp;
536 In.segmentlist[2*In_Index] = In_Index;
537 In.segmentlist[2*In_Index+1] = In_Index+1;
538 In.segmentmarkerlist[In_Index] = CurrComp;
544 In.segmentlist[2*(In_Index-1)+1] = 0;
555 phi1 = 0.000000E+0000;
560 s = (phi2- phi1)/(
double)N_Interface_Vert;
566 for(i=0;i<N_Interface_Vert;i++)
568 theta = phi1 + (double)i*s;
572 In.pointlist[2*In_Index] = T_a*cos(theta);;
573 In.pointlist[2*In_Index+1] = T_b*sin(theta);
589 In.pointmarkerlist[In_Index] = CurrComp;
590 In.segmentlist[2*In_Index] = In_Index;
591 In.segmentlist[2*In_Index+1] = In_Index+1;
592 In.segmentmarkerlist[In_Index] = CurrComp;
596 In.segmentlist[2*(In_Index-1)+1] = temp;
600 if(Out.pointlist!=NULL) {
601 free(Out.pointlist); Out.pointlist = NULL;}
602 if(Out.pointattributelist!=NULL) {
603 free(Out.pointattributelist); Out.pointattributelist = NULL;}
604 if(Out.pointmarkerlist!=NULL) {
605 free(Out.pointmarkerlist); Out.pointmarkerlist = NULL;}
606 if(Out.trianglelist!=NULL) {
607 free(Out.trianglelist); Out.trianglelist = NULL;}
608 if(Out.triangleattributelist!=NULL) {
609 free(Out.triangleattributelist); Out.triangleattributelist = NULL;}
610 if(Out.trianglearealist!=NULL) {
611 free(Out.trianglearealist); Out.trianglearealist = NULL;}
612 if(Out.neighborlist!=NULL) {
613 free(Out.neighborlist); Out.neighborlist = NULL;}
614 if(Out.segmentlist!=NULL) {
615 free(Out.segmentlist); Out.segmentlist = NULL;}
616 if(Out.segmentmarkerlist!=NULL) {
617 free(Out.segmentmarkerlist); Out.segmentmarkerlist = NULL;}
618 if(Out.holelist!=NULL) {
619 free(Out.holelist); Out.holelist = NULL;}
620 if(Out.regionlist!=NULL) {
621 free(Out.regionlist); Out.regionlist = NULL;}
622 if(Out.edgelist!=NULL) {
623 free(Out.edgelist); Out.edgelist = NULL;}
624 if(Out.edgemarkerlist!=NULL) {
625 free(Out.edgemarkerlist); Out.edgemarkerlist = NULL;}
626 if(Out.normlist!=NULL) {
627 free(Out.normlist); Out.normlist = NULL;}
630 triangulate((
char*)opts.str().c_str(), &In, &Out, (
struct triangulateio *)NULL);
638 VertexDel =
new TVertex*[3*N_RootCells];
642 for(i=0;i<N_Cells;i++)
647 for(j=0;j<N_Joints;j++)
651 VertexDel[CurrVertex] = cell->
GetVertex(j);
657 for(k=0;k<CurrVertex;k++)
664 VertexDel[CurrVertex] = cell->
GetVertex(j);
670 for(k=0;k<CurrVertex;k++)
671 if(VertexDel[k]==cell->
GetVertex((j+1)%N_Vertices))
677 VertexDel[CurrVertex] = cell->
GetVertex((j+1)%N_Vertices);
682 for(i=0;i<CurrVertex;i++)
686 OutPut(CurrVertex<<
" vertices were deleted"<<endl);
689 for(i=0;i<N_RootCells;i++)
691 OutPut(N_RootCells<<
" cells were deleted"<<endl);
698 UpdateBound[0]->
SetParams(Xi[0], Yi[0], Xi[1]-Xi[0],Yi[1]-Yi[0]);
699 UpdateBound[1]->
SetParams(Xi[1], Yi[1], Xi[2]-Xi[1],Yi[2]-Yi[1]);
700 UpdateBound[2]->
SetParams(Xi[2], Yi[2], Xi[3]-Xi[2],Yi[3]-Yi[2]);
701 UpdateBound[3]->
SetParams(Xi[3], Yi[3], Xi[0]-Xi[3],Yi[0]-Yi[3]);
705 UpdateIntface->
SetParams(C_x, C_y, T_a, T_b, phi1, phi2);
708 N_RootCells = Out.numberoftriangles;
711 Coordinates = Out.pointlist;
712 Triangles = Out.trianglelist;
713 PartMarker =
new int[Out.numberofpoints];
716 N_G = Out.numberofpoints;
717 NewVertices =
new TVertex*[N_G];
720 NewVertices[i] =
new TVertex(Coordinates[2*i], Coordinates[2*i+1]);
746 for (i=0;i<N_RootCells;i++)
750 CellTree[i]->
SetVertex(0, NewVertices[Out.trianglelist[3*i ]]);
751 CellTree[i]->
SetVertex(1, NewVertices[Out.trianglelist[3*i + 1]]);
752 CellTree[i]->
SetVertex(2, NewVertices[Out.trianglelist[3*i + 2]]);
754 ((
TMacroCell *) CellTree[i])->SetSubGridID(0);
767 N_G = Out.numberofpoints;
768 PointNeighb =
new int[N_G];
770 memset(PointNeighb, 0, N_G *SizeOfInt);
772 for (i=0;i<3*N_RootCells;i++)
773 PointNeighb[Triangles[i]]++;
778 if (PointNeighb[i] > maxEpV) maxEpV = PointNeighb[i];
779 delete [] PointNeighb;
781 PointNeighb =
new int[++maxEpV * N_G];
783 memset(PointNeighb, 0, maxEpV * N_G *SizeOfInt);
788 for(i=0;i<3*N_RootCells;i++)
790 j = Triangles[i]*maxEpV;
792 PointNeighb[j + PointNeighb[j]] = i / 3;
796 N_G = Out.numberofedges;
799 a = Out.edgelist[2*i];
800 b = Out.edgelist[2*i+1];
805 len1 = PointNeighb[a*maxEpV];
806 len2 = PointNeighb[b*maxEpV];
809 for (j=1;j<=len1;j++)
811 Neighb_tmp = PointNeighb[a*maxEpV + j];
812 for (k=1;k<=len2;k++)
813 if (Neighb_tmp == PointNeighb[b*maxEpV + k])
815 Neib[CurrNeib++] = Neighb_tmp;
818 if (CurrNeib == 2)
break;
825 if (Out.edgemarkerlist[i])
827 CurrComp = Out.edgemarkerlist[i] - 1;
828 if (CurrComp >= 100000) CurrComp -= 100000;
836 NewVertices[a]->GetX(), NewVertices[a]->GetY(), T_a) ||
838 NewVertices[b]->GetX(), NewVertices[b]->GetY(), T_b))
840 cerr<<
"Error: could not set parameter values"<<endl;
841 OutPut(NewVertices[a]<<endl);
842 OutPut(NewVertices[b]<<endl);
854 if(BDpart==1 && CurrComp==0 && fabs(T_a)==0 ) T_a=1;
865 T_a, T_b, CellTree[Neib[0]], CellTree[Neib[1]]);
870 T_a, T_b, CellTree[Neib[0]], CellTree[Neib[1]]);
889 cerr <<
"Error!!!!!!!! not enough neighbours!" << endl;
891 Joint =
new TJointEqN(CellTree[Neib[0]], CellTree[Neib[1]]);
896 if (Triangles[3*Neib[0]+j] == a)
break;
900 if (Triangles[3*Neib[0]+k] == b)
break;
920 CellTree[Neib[0]]->
SetJoint(j, Joint);
926 if (Triangles[3*Neib[1]+j] == a)
break;
930 if (Triangles[3*Neib[1]+k] == b)
break;
950 CellTree[Neib[1]]->
SetJoint(j, Joint);
953 if (Joint->
GetType() == InterfaceJoint ||
954 Joint->
GetType() == IsoInterfaceJoint)
959 delete [] NewVertices;
960 delete [] PointNeighb;
961 delete [] In.pointlist;
962 delete [] In.pointmarkerlist;
963 delete [] In.segmentlist;
964 delete [] In.segmentmarkerlist;
966 if(Out.pointlist!=NULL) {
967 free(Out.pointlist); Out.pointlist = NULL;}
968 if(Out.pointattributelist!=NULL) {
969 free(Out.pointattributelist); Out.pointattributelist = NULL;}
970 if(Out.pointmarkerlist!=NULL) {
971 free(Out.pointmarkerlist); Out.pointmarkerlist = NULL;}
972 if(Out.trianglelist!=NULL) {
973 free(Out.trianglelist); Out.trianglelist = NULL;}
974 if(Out.triangleattributelist!=NULL) {
975 free(Out.triangleattributelist); Out.triangleattributelist = NULL;}
976 if(Out.trianglearealist!=NULL) {
977 free(Out.trianglearealist); Out.trianglearealist = NULL;}
978 if(Out.neighborlist!=NULL) {
979 free(Out.neighborlist); Out.neighborlist = NULL;}
980 if(Out.segmentlist!=NULL) {
981 free(Out.segmentlist); Out.segmentlist = NULL;}
982 if(Out.segmentmarkerlist!=NULL) {
983 free(Out.segmentmarkerlist); Out.segmentmarkerlist = NULL;}
984 if(Out.holelist!=NULL) {
985 free(Out.holelist); Out.holelist = NULL;}
986 if(Out.regionlist!=NULL) {
987 free(Out.regionlist); Out.regionlist = NULL;}
988 if(Out.edgelist!=NULL) {
989 free(Out.edgelist); Out.edgelist = NULL;}
990 if(Out.edgemarkerlist!=NULL) {
991 free(Out.edgemarkerlist); Out.edgemarkerlist = NULL;}
992 if(Out.normlist!=NULL) {
993 free(Out.normlist); Out.normlist = NULL;}
double * GetValues()
Definition: FEFunction2D.h:67
virtual int GetTofXY(double X, double Y, double &T)=0
static TFE2D * GetFE2D(FE2D FE)
Definition: FEDatabase2D.h:353
int GetLength()
Definition: FEFunction2D.h:63
void GetTreeInfo(TBaseCell **&celltree, int &N_rootcells)
get tree of cells
Definition: Domain.h:176
TCollection * GetCollection() const
Definition: FESpace.h:131
JointType GetType()
Definition: Joint.h:75
TBaseCell * GetCell(int i) const
return Cell with index i in Cells-array
Definition: Collection.h:50
TBoundComp2D * GetBoundComp() const
Definition: BoundEdge.h:77
int GetBdPartID(int BdCompID)
get boundary part of BdCompID
Definition: Domain.C:88
void SetParams(double xmid, double ymid, double radius_a, double radius_b, double phi1, double phi2)
Definition: BdCircle.C:21
Definition: IsoBoundEdge.h:18
TBoundComp2D * GetBdComp(int i)
Definition: BoundPart.h:49
int SetJoint(int J_i, TJoint *J)
set the pointer to face J_i to J
Definition: BaseCell.h:168
contains the boundary description, the virtual cell tree and macro grid
Definition: Domain.h:36
void SetParams(double xstart, double ystart, double delx, double dely)
Definition: BdLine.C:21
Definition: FESpace2D.h:28
double RE_NR
Definition: Database.h:313
int SetParam(TDomain *domain)
Definition: Iterator.C:17
static double ** GetOrigElementValues(BaseFunct1D BaseFunct, MultiIndex1D MultiIndex)
Definition: FEDatabase2D.h:300
static TIterator ** IteratorDB
Definition: Database.h:1131
store cells in an array, used by cell iterators
Definition: Collection.h:18
int GetN_Joints()
return the number of joints
Definition: BaseCell.h:185
virtual TVertex * GetVertex(int Vert_i)=0
return the pointer to vertex with number i
represent a unit of the macro grid
Definition: MacroCell.h:15
Definition: JointEqN.h:20
int GetN_Cells() const
return number of cells
Definition: Collection.h:46
TJoint * GetJoint(int J_i)
return the pointer to face with number i
Definition: BaseCell.h:175
TFESpace2D * GetFESpace2D()
Definition: FEFunction2D.h:59
TBoundPart * GetBdPart(int i)
get i-th boundary part
Definition: Domain.h:172
FE2D GetFE2D(int i, TBaseCell *cell)
Definition: FESpace2D.C:1184
TFEDesc2D * GetFEDesc2D() const
Definition: FE2D.h:97
int GetN_Vertices()
return the number of vertices of the cell
Definition: BaseCell.h:179
int ** GetJointDOF() const
Definition: FEDesc2D.h:86
static BaseFunct2D * GetBaseFunct2D_IDFromFE2D()
Definition: FEDatabase2D.h:417
Definition: BoundPart.h:21
void SetTreeInfo(TBaseCell **celltree, int N_rootcells)
set tree of cells
Definition: Domain.h:183
static int * GetN_BaseFunctFromFE2D()
Definition: FEDatabase2D.h:421
bool IsFreeBoundary() const
Definition: BoundComp.h:57
TCollection * GetCollection(Iterators it, int level)
produce a collection with all cells returned by iterator it
Definition: Domain.C:1982
int GetID() const
Definition: BoundComp.h:49
Definition: BdCircle.h:18
information for finite element data structure
Definition: BaseCell.h:25
int * GetGlobalNumbers() const
Definition: FESpace.h:135
int GetN_Edges()
return the number of edges of the cell
Definition: BaseCell.h:182
static TRefDesc ** RefDescDB
Definition: Database.h:1125
Definition: BoundEdge.h:19
Definition: IsoInterfaceJoint.h:18
Definition: BoundComp.h:27
static RefTrans2D GetOrig(int N_LocalUsedElements, FE2D *LocalUsedElements, TCollection *Coll, TBaseCell *cell, bool *Needs2ndDer, int &N_Points, double *&xi, double *&eta, double *&weights, double *X, double *Y, double *absdetjk)
Definition: FEDatabase2D.C:1765
virtual int SetVertex(int Vert_i, TVertex *Vert)=0
set the pointer of vertex Vert_i to Vert
Definition: FEDesc2D.h:15
represent geometric information of the cell
Definition: GridCell.h:15
int * GetBeginIndex() const
Definition: FESpace.h:142
int GetN_JointDOF() const
Definition: FEDesc2D.h:61
Definition: InterfaceJoint.h:18
static TParamDB * ParamDB
Definition: Database.h:1134
int GetLocalBdCompID(int BdCompID)
get local number of boundary component
Definition: Domain.C:98
Definition: FEFunction2D.h:24