00001 #include "meshtransfert.h" 00002 #include "globalt.h" 00003 #include "globmatt.h" 00004 #include "gtopology.h" 00005 #include "gelement.h" 00006 #include "gnode.h" 00007 #include "gmatrix.h" 00008 #include "mathem.h" 00009 #include <time.h> 00010 #include <math.h> 00011 #include <stdlib.h> 00012 #include "least_square.h" 00013 /** 00014 Function approximates values stored on integration points to nodes. 00015 Used interpolation method is 'spr_smoothing'. 00016 Values in nodes are returned by array nodvalues. 00017 Value position in array is (val on nod[0] ... val on nod[gt->nn]) 00018 gt->nadjelem have to be allocated == to use function gt->adjacelem. 00019 00020 @param ipvalues - values in integration points - array, dimension is total number of ip (tinp) 00021 00022 @param nodvalues - empty(returned) array, dimension is number of nodes (Tt->nn) 00023 00024 created 9.12.2002, Ladislav Svoboda, termit@cml.fsv.cvut.cz 00025 modified by Tomas Krejci, 13.11.2003 00026 */ 00027 /* 00028 void give_nodvals_ip (double *ipvalues,double *nodvalues) 00029 { 00030 long i,j,k,nip,ipp,ntm,*nsp,dim; 00031 double **spcoord,**spvalue,**ipcoord; 00032 00033 nsp = new long [Tt->ne]; 00034 spcoord = new double* [Tt->ne]; 00035 spvalue = new double* [Tt->ne]; 00036 00037 ntm = Tp->ntm; 00038 00039 ipcoord = new double* [Tm->tnip]; 00040 for (i=0;i<Tm->tnip;i++) 00041 ipcoord[i] = new double [3]; 00042 00043 allipcoord (ipcoord); 00044 00045 for (i=0;i<Tt->ne;i++){ 00046 nip = 0; 00047 for (j=0;j<ntm;j++) 00048 for (k=0;k<ntm;k++) 00049 nip += Tt->give_nip(i,j,k); 00050 00051 nsp[i] = nip; 00052 dim = Tt->give_dimension (i); 00053 00054 ipp = Tt->elements[i].ipp[0][0]; 00055 spcoord[i] = new double [nip*dim]; 00056 spvalue[i] = new double [nip]; 00057 for (j=0;j<nip;j++){ 00058 for (k=0;k<dim;k++) 00059 spcoord[i][j*dim+k] = ipcoord[ipp+j][k]; 00060 spvalue[i][j] = ipvalues[ipp+j]; 00061 } 00062 } 00063 00064 least_square spr_est(dim,1,Gtt,1); 00065 00066 spr_est.run (Outt,nsp,spcoord,spvalue,nodvalues); 00067 00068 delete [] nsp; 00069 for (i=0;i<Tt->ne;i++){ delete [] spcoord[i]; delete [] spvalue[i]; } 00070 delete [] spcoord; 00071 delete [] spvalue; 00072 for (i=0;i<Tm->tnip;i++) delete [] ipcoord[i]; 00073 delete [] ipcoord; 00074 } 00075 */ 00076 00077 /** 00078 function computes coordinates of all integration points in "global" mesh 00079 00080 @param ipcoord - empty array, dimension is Tt->tnip x 3 00081 00082 created 5.12.2002, Ladislav Svoboda, termit@cml.fsv.cvut.cz 00083 modified by Tomas Krejci 13.11.2003 00084 */ 00085 /* 00086 void allipcoord (double **ipcoord) 00087 { 00088 long i,ri,ci,ntm,ipp; 00089 00090 for (i=0;i<Tt->ne;i++){ 00091 ntm = Tp->ntm; 00092 for (ri=0;ri<ntm;ri++){ 00093 for (ci=0;ci<ntm;ci++){ 00094 ipp = Tt->elements[i].ipp[ri][ci]; 00095 00096 switch (Tt->give_elem_type(i)){ 00097 //case barlint:{ Lbt->ipcoordblock (i,ri,ci,ipcoord+ipp); break; } 00098 //case barlintax:{ Lbat->ipcoordblock (i,ri,ci,ipcoord+ipp); break; } 00099 //case barquadt:{ Qbt->ipcoordblock (i,ri,ci,ipcoord+ipp); break; } 00100 //case barquadtax:{ Qbat->ipcoordblock (i,ri,ci,ipcoord+ipp); break; } 00101 case trlint:{ Ltt->ipcoordblock (i,ri,ci,ipcoord+ipp); break; } 00102 case trlaxisym:{ Ltat->ipcoordblock (i,ri,ci,ipcoord+ipp); break; } 00103 case quadlint:{ Lqt->ipcoordblock (i,ri,ci,ipcoord+ipp); break; } 00104 case quadquadt:{ Qqt->ipcoordblock (i,ri,ci,ipcoord+ipp); break; } 00105 case quadquadtax:{ Qqat->ipcoordblock (i,ri,ci,ipcoord+ipp); break; } 00106 case quadlaxisym:{ Lqat->ipcoordblock (i,ri,ci,ipcoord+ipp); break; } 00107 default:{ 00108 fprintf (stderr,"\n\n unknown element type is required in function allipcoord (%s, line %d).\n",__FILE__,__LINE__); 00109 } 00110 } 00111 } 00112 } 00113 } 00114 } 00115 */ 00116 00117 00118 /** 00119 function finds out node(from gt), which is the closest by required point 00120 00121 @param x,y,z - coordinates of point 00122 00123 created 30.1.2003, Ladislav Svoboda, termit@cml.fsv.cvut.cz 00124 */ 00125 /* 00126 long adjacnodet (gtopology *gt,double x,double y,double z) 00127 { 00128 long i,a = 0; 00129 double d,min = 1e50; 00130 00131 for(i=0;i<gt->ne;i++){ 00132 d = (gt->gnodes[i].x - x)*(gt->gnodes[i].x - x) + (gt->gnodes[i].y - y)*(gt->gnodes[i].y - y) + (gt->gnodes[i].z - z)*(gt->gnodes[i].z - z); 00133 if (d<min){ 00134 min=d; 00135 a=i; 00136 } 00137 } 00138 return (a); 00139 } 00140 00141 */