ELEMENTO TRIANGULAR DE TRANSICION - 4 NODOS

Implementación en Mathematica

v.2018

1. DATOS INICIALES

INICIO

Off[General::"spell1"]
Off[General::"spell"]

SetDirectory[NotebookDirectory[]]

DEFINICION  ELEMENTO TRIANGULAR DE TRANSICION DE 4 NODOS

DEFINICION GRAFICA

TriT4=;

TriT4r=Show[TriT4, ImageSize->300]

COORDENADAS TRIANGULARES NODOS - ###

Cnt={{1,0,0},{0,1,0},X,{(1/2),(1/2),0}};

NNodos=Dimensions[Cnt][[1]]

ELEMENTO COMPLETO NECESARIO - 2 DIVISIONES POR LADO

TriR6=;

TriR6r=Show[TriR6, ImageSize->300]

DEFINICION  ELEMENTO BASE REAL - COORDENADAS CARTESIANAS

COORDENADAS CARTESIANAS NODOS ELEMENTO BASE REAL

NNodosB=3;

Cne=Table[{0,0},{i,NNodosB}];

Cne[[1]]={-10,0};Cne[[2]]={10,0};Cne[[3]]={0,10*Sqrt[3]};

FUNCION TRANSFORMACION DE COORDENADAS DE TRIANGULARES A CARTESIANAS

TtC[ζ1_,ζ2_,ζ3_]=(({{1, 1, 1}, {Cne[[1]][[1]], Cne[[2]][[1]], Cne[[3]][[1]]}, {Cne[[1]][[2]], Cne[[2]][[2]], Cne[[3]][[2]]}})).(({{ζ1}, {ζ2}, {ζ3}}));

COORDENADAS CARTESIANAS DE LOS NODOS

Cnc=Table[{0,0},{i,NNodos}];

Do[
Cnc[[i]]={TtC[Cnt[[i]][[1]],Cnt[[i]][[2]],Cnt[[i]][[3]]][[2,1]],TtC[Cnt[[i]][[1]],Cnt[[i]][[2]],Cnt[[i]][[3]]][[3,1]]},
{i,NNodos}
];

IMAGEN DEL ELEMENTO - #

FUNCION REPRESENTACION GRAFICA ELEMENTO Y NODOS

ElementPlot[b_List,options___]:=Module[{asa,color,nr,circles,lines},asa=Select[{options},((!SameQ[#[[1]],NodeColor])&&(!SameQ[#[[1]],NodeSize]))&];
{color,nr}={NodeColor,NodeSize}/.{options}/.{NodeColor->GrayLevel[0],NodeSize->PointSize[0.06]};
circles=Map[Point[#]&,Partition[Flatten[b],2]];
lines=Line[Append[b[[1]],First[b[[1]]]]];
Show[Graphics[{nr,color,circles}],Graphics[lines],Evaluate[Sequence[##]&@@asa]]];

DEFINICION VECTOR DE NODOS- ###

ptsexteriores={Cnc[[1]],Cnc[[4]],X,Cnc[[3]]};

ptsinteriores={};

IMAGEN DE COMPROBACION

Imagen=Elemento=ElementPlot[{ptsexteriores,ptsinteriores},AspectRatio->Automatic,Axes->True,PlotRange->{{-12,12},{-2,20}},ImageSize-> 250,NodeColor->RGBColor[0,0,1]]

3. FUNCIONES DE FORMA - METODO PRODUCTO DE CURVAS - #

CURVAS A CONSIDERAR - TRIANGULO REGULAR DE 6 NODOS

TriR6r

Cu=Table[0,{i,9}];

LADOS

Cu[[1]]=ζ3;Cu[[2]]=ζ1;Cu[[3]]=ζ2;

MEDIANAS

Cu[[4]]=(ζ1-1/2);

Cu[[5]]=(ζ2-1/2);

Cu[[6]]=(ζ3-1/2);

DEFINICION PRODUCTO DE CURVAS EN CADA NODO - NODOS NO ESQUINA

Nc=Table[0,{i,NNodos}];

Show[TriT4,ImageSize-> 350]

Tipo 2 - LADOS

Nc[[4]]=Cu[[2]]*Cu[[3]];

OBTENCION FUNCIONES DE FORMA - NODOS NO ESQUINA

Clear[Nf]

Nfp=Table[0,{i,NNodos}];

Nf=Table[0,{i,NNodos}];

Do[
Nfp[[i]]=a*Nc[[i]];
eq=1==Nfp[[i]]/.{ζ1 ->Cnt[[i,1]],ζ2->Cnt[[i,2]],ζ3->Cnt[[i,3]]};
as=a/.Solve[eq,a][[1]];Print["Nodo ",i];
Nf[[i]]=Simplify[Nfp[[i]]/.{a->as}],
{i,4,NNodos}
];

MatrixForm[Nf]

OBTENCION FUNCIONES DE FORMA - NODOS ESQUINA

Utilizamos las Funciones de Forma del Triangulo de 3 Nodos.

Nf3={ζ1,ζ2,ζ3};

NODO 1

Clear[a4]

Nf[[1]]=Nf3[[1]]+a4*Nf[[4]];

Show[TriT4,ImageSize-> 400]

eq=0==Nf[[1]]/.{ζ1->Cnt[[4]][[1]],ζ2->Cnt[[4]][[2]],ζ3->Cnt[[4]][[3]]}
a4s=a4/.Solve[eq,a4][[1]]

Nf[[1]]=Simplify[Nf[[1]]/.{a4->a4s}]

NODO 2

Clear[a4]

Nf[[2]]=Nf3[[2]]+a4*Nf[[4]];

Show[TriT4,ImageSize-> 400]

eq=0==Nf[[2]]/.{ζ1->Cnt[[4]][[1]],ζ2->Cnt[[4]][[2]],ζ3->Cnt[[4]][[3]]}
a4s=a4/.Solve[eq,a4][[1]]

Nf[[2]]=Simplify[Nf[[2]]/.{a4->a4s}]

NODO 3

Clear[a4]

Nf[[3]]=Nf3[[3]]+a4*Nf[[4]];

Show[TriT4,ImageSize-> 400]

eq=0==Nf[[3]]/.{ζ1->Cnt[[4]][[1]],ζ2->Cnt[[4]][[2]],ζ3->Cnt[[4]][[3]]}
a4s=a4/.Solve[eq,a4][[1]]

Nf[[3]]=Simplify[Nf[[3]]/.{a4->a4s}]

FUNCIONES DE FORMA DE TODOS LOS NODOS

MatrixForm[Nf]

COMPROBACION SUMA UNIDAD

Suma=(∑_(i=1)^NNodos)Nf[[i]]

Simplify[Suma/.{ζ1->1-ζ2-ζ3}]

OK - SE CUMPLE LA CONDICION DE COMPLETITUD

Proceso para comprobar Valor Funciones de Forman en Nodos - en caso de Error

(*Do[
Print["NODO ",j];
Do[
Print[i,"  ",Simplify[Nf[[j]]/.{ζ1 ->Cnt[[i,1]],ζ2->Cnt[[i,2]],ζ3->Cnt[[i,3]]}]/.{ζ1->1-ζ2-ζ3}],
{i,NNodos}
  ],
{j,NNodos}
];*)

REPRESENTACION GRAFICA

Función Representación Gráfica Funciones de Forma

PlotTriangleShapeFunction[xytrig_,f_,Nsub_,aspect_]:=Module[{Ni,line3D={},poly3D={},zc1,zc2,zc3,xyf1,xyf2,xyf3,xc,yc,x1,x2,x3,y1,y2,y3,z1,z2,z3,iz1,iz2,iz3,d},{{x1,y1,z1},{x2,y2,z2},{x3,y3,z3}}=Take[xytrig,3];xc={x1,x2,x3};yc={y1,y2,y3};Ni=Nsub*3;Do[Do[iz3=Ni-iz1-iz2;If[iz3≤0,Continue[]];d=0;If[Mod[iz1+2,3]==0&&Mod[iz2-1,3]==0,d=1];If[Mod[iz1-2,3]==0&&Mod[iz2+1,3]==0,d=-1];If[d==0,Continue[]];zc1=N[{iz1+d+d,iz2-d,iz3-d}/Ni];zc2=N[{iz1-d,iz2+d+d,iz3-d}/Ni];zc3=N[{iz1-d,iz2-d,iz3+d+d}/Ni];xyf1={xc.zc1,yc.zc1,f[zc1[[1]],zc1[[2]],zc1[[3]]]};xyf2={xc.zc2,yc.zc2,f[zc2[[1]],zc2[[2]],zc2[[3]]]};xyf3={xc.zc3,yc.zc3,f[zc3[[1]],zc3[[2]],zc3[[3]]]};AppendTo[poly3D,Polygon[{xyf1,xyf2,xyf3}]];AppendTo[line3D,Line[{xyf1,xyf2,xyf3,xyf1}]],{iz2,1,Ni-iz1}],{iz1,1,Ni}];Show[Graphics3D[RGBColor[1,0,0]],Graphics3D[poly3D],Graphics3D[Thickness[.002]],Graphics3D[line3D],Graphics3D[RGBColor[0,0,0]],Graphics3D[Thickness[.005]],Graphics3D[Line[xytrig]],PlotRange->All,BoxRatios->{1,1,aspect},Boxed->False]];

Representación Gráfica Funciones Forma Elemento.

Ng=Table[0,{i,NNodos}];

xyc1={0,0,0};xyc2={3,0,0};xyc3={Sqrt[3],3/2,0};xytrig=N[{xyc1,xyc2,xyc3,xyc1}];

Control de Cuadricula

Nsub=15;

Do[
fi[ζ1_,ζ2_,ζ3_]=Nf[[i]];
Ng[[i]]=PlotTriangleShapeFunction[xytrig,fi,Nsub,1/2],
{i,NNodos}
];

4. RESULTADOS INTERACTIVOS - #

Manipulate[{TriT4r,Ng[[n]],Nf[[n]]},{n,1,Dimensions[Nf][[1]],1},{n,Range[Dimensions[Nf][[1]]]},FrameLabel->{"FUNCION DE FORMA EN NODO n - TRIANGULO TRANSICION 4 NODOS"}, SaveDefinitions->True]

5. DERIVADAS FUNCIONES DE FORMA Y JACOBIANO - #

Inicializaciones Necesarias - #

ClearAll[x1,x2,x3,x4,y1,y2,y3,y4,ζ1,ζ2,ζ3];

1 - Definción Isoparamétrica del Elemento - ###

IsoP=;

IsoPr=Show[IsoP,ImageSize-> 650]

x={x1,x2,x3,X};
y={y1,X,y3,y4};

2 - Funciones de Forma  <<< -----

Nf

3 - Derivadas Funciones de Forma respecto Coordenadas Naturales - ###

Nf1=D[Nf,ζ1];
Nf2=D[Nf,ζ2];
Nf3=X;
{Nf1,Nf2,Nf3}=Simplify[{Nf1,Nf2,Nf3}];

Nf1//MatrixForm

Nf2//MatrixForm

Nf3//MatrixForm

4 - Derivadas Coordenadas Triangulares respecto a las Cartesianas - Desarrollo

Calculo Elementos Matriz Jacobiana - Elemento Considerado - ###

SistemA=;

SistemAr=Show[SistemA,ImageSize-> 650]

{Nf1x,Nf2x,Nf3x}=Simplify[{x.Nf1,x.Nf2,X}];
{X,Nf2y,Nf3y}=Simplify[{y.Nf1,y.Nf2,y.Nf3}];

{Nf1x,Nf2x,Nf3x} //MatrixForm

{Nf1y,Nf2y,Nf3y} //MatrixForm

JacobianA=;

JacobianAr=Show[JacobianA,ImageSize-> 650]

Para mostrar el planteamiento genérico, de momento no aplicamos estas definiciones:

Clear[Jx1,Jx2,Jx3,Jy1,Jy2,Jy3]

(*Jx1=Nf1x; Jx2=Nf2x; Jx3=Nf3x;*)

(*Jy1=Nf1y; Jy2=Nf2y; Jy3=Nf3y;*)

Definición Matriz Jacobiana y Sistema de Ecuaciones a Resolver - Planteamiento Genérico - ###

J=(({{1, 1, 1}, {Jx1, X, Jx3}, {Jy1, Jy2, Jy3}}));

Determinante de la Matriz Jacobiana - J

Jdet=Det[J]

Defincin de Jc

Jc = 1/2*Jdet

Matriz Derivadas Coordenadas Triangulares respecto a las Cartesianas - Matriz Incognitas - ###

P=(({{f1x, f1y}, {f2x, f2y}, {f3x, X}}));

Matriz de Terminos Independientes - ###

Ti=(({{X, 0}, {1, 0}, {0, 1}}));

Solución del Sistema de Ecuaciones

P=Inverse[J].Ti

Definción Matriz Jacobiana Inversa Modificada - Jim - ###

JacobianaModificadA=;

JacobianaModificadAr=Show[JacobianaModificadA,ImageSize-> 650]

Jc

P

Jim=2*Jc*P

Jim // MatrixForm

(({{Jy23, X}, {Jy31, Jx13}, {Jy12, Jx21}}))=Jim;

Jx13

Necesario para mantener el planteamiento genérico:

Clear[Jdet,Jc,Jy23,Jx32,Jy31,Jx13,Jy12,Jx21]

Jc = 1/2*Jdet;

JacobianaModificadAr

(({{f1x, f1y}, {f2x, X}, {f3x, f3y}}))=(1/(2*Jc))*(({{Jy23, Jx32}, {Jy31, Jx13}, {X, Jx21}}));

f1x

5 - Derivadas de las Funciones de Forma respecto a las Coordenadas Cartesianas - Desarrollo

Definición de las Derivadas de las Funciones de Forma respecto a las Coordenadas Cartesianas - Planteamiento Genérico - ###

DerivadasFuncionesFormaCartesianaS=;

DerivadasFuncionesFormaCartesianaSr=Show[DerivadasFuncionesFormaCartesianaS,ImageSize-> 650]

DerivadasFuncionesFormaCartesianasM=;

DerivadasFuncionesFormaCartesianasMr=Show[DerivadasFuncionesFormaCartesianasM,ImageSize-> 650]

(({{X}, {dNy}}))=Transpose[(({{f1x, X}, {f2x, f2y}, {f3x, f3y}}))].(({{Nf1}, {Nf2}, {X}}));

dNx & dNy  <<< -----

dNx

dNy

Se observa quedan en función de los elementos de la Matriz Jacobiana Inversa Modificada - Jim

Matriz Jacobiana Inversa Modificada - Jim  - Planteamiento Genérico - ###

JacobianaModificadAr

(({{Jy23, X}, {Jy31, Jx13}, {Jy12, Jx21}}))=Jim;

Jx13

Elemetos Matriz Jacobiana - Elemento Considerado - ###

JacobianAr

Aplicamos las deiniciones correspondientes al alemento considerado:

Jx1=Nf1x; Jx2=Nf2x; Jx3=Nf3x;

Jy1=Nf1y; Jy2=X; Jy3=Nf3y;

Matriz Jacobiana Inversa Modificada - Jim - Elemento Considerado

(({{Jy23, X}, {Jy31, Jx13}, {Jy12, Jx21}}));

Jy23

Jx32

Jy31

Jx13

Jy12

Jx21

MODULO CALCULO FUNCIONES DE FORMA Y JACOBIANO

MODULO GENERICO A COMPLETAR - #

MODULO GENERICO A COMPLETAR: XX = No. Nodos, {xI,yI}= nodo i-esimo

(*TrigXXIsoPShapeFunDer[ncoor_,tcoor_]:=Module[{ζ1,ζ2,ζ3,x1,x2,x3,x4,x5,x6,x7,x8,x9,xI,y1,y2,y3,y4,y5,y6,y7,y8,y9,yI,Jx21,Jx32,Jx13,Jy12,Jy23,Jy31,Nf,dNx,dNy,Jdet},{ζ1,ζ2,ζ3}=tcoor;
{{x1,y1},{x2,y2},{x3,y3},{x4,y4},{x5,y5},{x6,y6},{x7,y7},{x8,y8},{x9,y9},{xI,yI}}=ncoor;
Nf=(*Nf*);
Jx21=(*Jx21*);
Jy12=(*Jy12*);
Jx13=(*Jx13*);
Jy31=(*Jy31*);
Jx32=(*Jx32*);
Jy23=(*Jy23*);
Jdet=Jx21*Jy31-Jy12*Jx13;
dNx=(*dNx*);
dNy=(*dNy*);
Return[Simplify[{Nf,dNx,dNy,Jdet}]]];*)

MODULO COMPLETADO - #

X;

6. MATRIZ DE RIGIDEZ - INTEGRACION NUMERICA - #

MODULO DE CALCULO DE LA MATRIZ DE RIGIDEZ

MODULO GENERICO A COMPLETAR - #

MODULO GENERICO: XX = No. Nodos, YY = Grados Libertad

(*TrigXXIsoPMembraneStiffness[ncoor_,mprop_,fprop_,opt_]:=Module[{i,k,l,p=3,numer=False,Emat,th={fprop},h,tcoor,w,c,Nf,dNx,dNy,Jdet,B,Ke=Table[0,{YY},{YY}]},Emat=mprop[[1]];If[Length[fprop]>0,th=fprop[[1]]];
If[Length[opt]>0,numer=opt[[1]]];
If[Length[opt]>1,p=opt[[2]]];
If[p≠1&&p≠-3&&p≠3&&p≠6&&p≠-6&&p≠7&&p≠12,Print["Illegal p"];Return[Null]];
For[k=1,k≤Abs[p],k++,{tcoor,w}=TrigGaussRuleInfo[{p,numer},k];
{Nf,dNx,dNy,Jdet}=TrigXXIsoPShapeFunDer[ncoor,tcoor];
If[Length[th]==0,h=th,h=th.Nf];c=w*Jdet*h/2;
B={Flatten[Table[{dNx[[i]],0},{i,XX}]],Flatten[Table[{0,dNy[[i]]},{i,XX}]],Flatten[Table[{dNy[[i]],dNx[[i]]},{i,XX}]]};
Ke+=c*Transpose[B].(Emat.B);];If[!numer,Ke=Simplify[Ke]];Return[Ke]];*)

MODULO COMPLETADO - #

X;

MODULO REGLAS DE CUADRATURA DE GAUSS

OPCION UNICA: DEFINICION DE CARLOS FELIPPA

TrigGaussRuleInfo[{rule_,numer_},point_]:=Module[{zeta,p=rule,i=point,g1,g2,g3,g4,w1,w2,w3,eps=10.^(-24),jkl={{1,2,3},{2,1,3},{1,3,2},{3,1,2},{2,3,1},{3,2,1}},info={{Null,Null,Null},0}},If[p==1,info={{1/3,1/3,1/3},1}];If[p==3,info={{1,1,1}/6,1/3};info[[1,i]]=2/3];If[p==-3,info={{1,1,1}/2,1/3};info[[1,i]]=0];If[p==6,g1=(8-Sqrt[10]+Sqrt[38-44*Sqrt[2/5]])/18;g2=(8-Sqrt[10]-Sqrt[38-44*Sqrt[2/5]])/18;If[i<4,info={{g1,g1,g1},(620+Sqrt[213125-53320*Sqrt[10]])/3720};info[[1,i]]=1-2*g1];If[i>3,info={{g2,g2,g2},(620-Sqrt[213125-53320*Sqrt[10]])/3720};info[[1,i-3]]=1-2*g2]];If[p==-6,If[i<4,info={{1,1,1}/6,3/10};info[[1,i]]=2/3];If[i>3,info={{1,1,1}/2,1/30};info[[1,i-3]]=0]];If[p==7,g1=(6-Sqrt[15])/21;g2=(6+Sqrt[15])/21;If[i<4,info={{g1,g1,g1},(155-Sqrt[15])/1200};info[[1,i]]=1-2*g1];If[i>3&&i<7,info={{g2,g2,g2},(155+Sqrt[15])/1200};info[[1,i-3]]=1-2*g2];If[i==7,info={{1/3,1/3,1/3},9/40}]];If[p==12,g1=0.063089014491502228340331602870819157;g2=0.249286745170910421291638553107019076;g3=0.053145049844816947353249671631398147;g4=0.310352451033784405416607733956552153;If[!numer,{g1,g2,g3,g4}=Rationalize[{g1,g2,g3,g4},eps]];w1=(30*g2^3*(4*g3^2+(1-2*g4)^2+4*g3*(-1+g4))+g3^2*(1-15*g4)+(-1+g4)*g4-g3*(-1+g4)*(-1+15*g4)+2*g2*(1+60*g3*g4*(-1+g3+g4))-6*g2^2*(3+10*(-1+g4)*g4+10*g3^2*(1+3*g4)+10*g3*(-1+g4)*(1+3*g4)))/(180*(g1-g2)*(-(g2*(-1+2*g2)*(-1+g3)*g3)+(-1+g3)*(g2-2*g2^2-2*g3+3*g2*g3)*g4-(g2*(-1+2*g2-3*g3)+2*g3)*g4^2+2*g1^2*(g2*(-2+3*g2)+g3-g3^2+g4-g3*g4-g4^2)+g1*(-4*g2^2+(-1+g3)*g3+(-1+g3)*(1+3*g3)*g4+(1+3*g3)*g4^2-2*g2*(-1+g3^2+g3*(-1+g4)+(-1+g4)*g4))));w2=(-1+12*(2-3*g1)*g1*w1+4*g3^2*(-1+3*w1)+4*g3*(-1+g4)*(-1+3*w1)+4*(-1+g4)*g4*(-1+3*w1))/(12*(g2*(-2+3*g2)+g3-g3^2+g4-g3*g4-g4^2));w3=(1-3*w1-3*w2)/6;If[i<4,info={{g1,g1,g1},w1};info[[1,i]]=1-2*g1];If[i>3&&i<7,info={{g2,g2,g2},w2};info[[1,i-3]]=1-2*g2];If[i>6,{j,k,l}=jkl[[i-6]];info={{0,0,0},w3};info[[1,j]]=g3;info[[1,k]]=g4;info[[1,l]]=1-g3-g4]];If[numer,Return[N[info]],Return[Simplify[info]]];];

7. TEST DEL ELEMENTO DE LADOS RECTOS SUPERPARAMETRICO.

DEFINICION DE LA GEOMETRIA

LadosRectosTesT=;

LadosRectosTesTr=Show[LadosRectosTesT,ImageSize-> 300]

DEFINICION DE LOS NODOS ELEMENTO - ###

ncoor=Cnc;

ptsexteriores={Cnc[[1]],Cnc[[4]],Cnc[[2]],X};

ptsinteriores={};

Elemento=ElementPlot[{ptsexteriores,ptsinteriores},AspectRatio->Automatic,Axes->True,PlotRange->{{-12,12},{-2,20}},NodeColor->RGBColor[0,0,1]]

DEFINICION DEL MATERIAL

ClearAll[Em,nu,h];

h=1;Em=288;nu=1/3;

Emat=Em/(1-nu^2)*{{1,nu,0},{nu,1,0},{0,0,(1-nu)/2}};

Print["Emat=",Emat//MatrixForm]

VERIFICACION DE LA MATRIZ DE RIGIDEZ - #

NUMERO DE PUNTOS DE GAUSS MINIMO PARA CONSEGUIR SUFICIENCIA DE RANGO - #

NF=NNodos*2.;

NG=((NF-3)/3)

Se necesitan como mínimo 2 Puntos -- Regla 3 minima

BUCLE GENERICO A COMPLETAR - #

BUCLE GENERICO: XX = No. Nodos, ZZ = GRADOS DE LIBERTAD - 3

(*For [i=1, i≤7,i++,p={1,-3,3,6,-6,7,12}[[i]];
Ke=TrigXXIsoPMembraneStiffness[ncoor,{Emat,0,0},{h},{True,p}];
Print["Gauss integration rule: ",p];
Print["Ke=",Chop[Simplify[Ke]]//MatrixForm];
Valores=Chop[Eigenvalues[N[Ke]]];
If[Valores[[ZZ]]!=0,Print["Valores propios matriz Ke=",Valores];Print["TENEMOS LA SUFICIENCIA DE RANGO PARA p=",p];Break[],Print["Valores proprios matriz Ke=",Valores];Print["NO tenemos la suficiencia de rango para p=",p]]
];*)

DESARROLLO DE LA MATRIZ DE RIGIDEZ - #

For [i=1, i≤7,i++,p={1,-3,3,6,-6,7,12}[[i]];
Ke=X[ncoor,{Emat,0,0},{h},{True,p}];
Print["Gauss integration rule: ",p];
Print["Ke=",Chop[Simplify[Ke]]//MatrixForm];
Valores=Chop[Eigenvalues[N[Ke]]];
If[Valores[[X]]!=0,Print["Valores propios matriz Ke=",Valores];Print["TENEMOS LA SUFICIENCIA DE RANGO PARA p=",p](*;Break[]*),Print["Valores proprios matriz Ke=",Valores];Print["NO tenemos la suficiencia de rango para p=",p]]
];

8. TEST DEL ELEMENTO DE LADOS CURVOS.

DEFINICION DE LA GEOMETRIA

LadosRectosTesT=;

LadosRectosTesTr=Show[LadosRectosTesT,ImageSize-> 300]

DEFINICION DE LOS NODOS ELEMENTO - #

ELEMENTO BASE REAL

Cng=Table[{0,0},{i,6}];

Cng[[1]]={-10,0};Cng[[2]]={10,0};Cng[[3]]={0,10*Sqrt[3]};Cng[[4]]={0,-10*Sqrt[3]*1/3};Cng[[5]]={10,10*Sqrt[3]*2/3};Cng[[6]]={-10,10*Sqrt[3]*2/3};

ptsexteriores={Cng[[1]],Cng[[4]],Cnc[[2]],Cng[[5]],Cng[[3]],Cng[[6]]};

ptsinteriores={};

Elemento=ElementPlot[{ptsexteriores,ptsinteriores},AspectRatio->Automatic,Axes->True,PlotRange->{{-12,12},{-10,20}},NodeColor->RGBColor[0,0,1]]

TRANSFORMACION DE COORDENADAS TRIANGULARES A CARTESIANAS

Nf6={ζ1*(2*ζ1-1),ζ2*(2*ζ2-1),ζ3*(2*ζ3-1),4*ζ1*ζ2,4*ζ2*ζ3,4*ζ3*ζ1};

TtC[ζ1_,ζ2_,ζ3_]=(({{1, 1, 1, 1, 1, 1}, {Cng[[1]][[1]], Cng[[2]][[1]], Cng[[3]][[1]], Cng[[4]][[1]], Cng[[5]][[1]], Cng[[6]][[1]]}, {Cng[[1]][[2]], Cng[[2]][[2]], Cng[[3]][[2]], Cng[[4]][[2]], Cng[[5]][[2]], Cng[[6]][[2]]}})).Nf6;

COORDENADAS TRIANGULARES DEL ELEMENTO CONSIDERADO

Cnt

COORDENADAS CARTESIANAS NODOS ELEMENTO REAL CONSIDERADO Y COMPROBACION GRAFICA - ###

Cnc=Table[{0,0},{i,NNodos}];

Do[
Cnc[[i]]={TtC[Cnt[[i]][[1]],Cnt[[i]][[2]],Cnt[[i]][[3]]][[2]],TtC[Cnt[[i]][[1]],Cnt[[i]][[2]],Cnt[[i]][[3]]][[3]]},
{i,NNodos}
];

Cnc

ncoor=Cnc;

ptsexteriores={Cnc[[1]],X,Cnc[[2]],Cnc[[3]]};

ptsinteriores={};

Elemento=ElementPlot[{ptsexteriores,ptsinteriores},AspectRatio->Automatic,Axes->True,PlotRange->{{-12,12},{-10,20}},NodeColor->RGBColor[0,0,1]]

DEFINICION DEL MATERIAL

ClearAll[Em,nu,h];

h=1;Em=7*72;nu=0;h=1;

Emat=Em/(1-nu^2)*{{1,nu,0},{nu,1,0},{0,0,(1-nu)/2}};

VERIFICACION DE LA MATRIZ DE RIGIDEZ - #

NUMERO DE PUNTOS DE GAUSS MINIMO PARA CONSEGUIR SUFICIENCIA DE RANGO - #

NF=NNodos*2.;

NG=((NF-3)/3)

Se necesitan como mínimo 2 Puntos -- Regla 3 minima

BUCLE GENERICO A COMPLETAR - #

BUCLE GENERICO: XX = No. Nodos, ZZ = GRADOS DE LIBERTAD - 3

(*For [i=1, i≤7,i++,p={1,-3,3,6,-6,7,12}[[i]];
Ke=TrigXXIsoPMembraneStiffness[ncoor,{Emat,0,0},{h},{True,p}];
Print["Gauss integration rule: ",p];
Print["Ke=",Chop[Simplify[Ke]]//MatrixForm];
Valores=Chop[Eigenvalues[N[Ke]]];
If[Valores[[ZZ]]!=0,Print["Valores propios matriz Ke=",Valores];Print["TENEMOS LA SUFICIENCIA DE RANGO PARA p=",p];Break[],Print["Valores proprios matriz Ke=",Valores];Print["NO tenemos la suficiencia de rango para p=",p]]
];*)

DESARROLLO DE LA MATRIZ DE RIGIDEZ - #

For [i=1, i≤7,i++,p={1,-3,3,6,-6,7,12}[[i]];
Ke=X[ncoor,{Emat,0,0},{h},{True,p}];
Print["Gauss integration rule: ",p];
Print["Ke=",Chop[Simplify[Ke]]//MatrixForm];
Valores=Chop[Eigenvalues[N[Ke]]];
If[Valores[[5]]!=0,Print["Valores propios matriz Ke=",Valores];Print["TENEMOS LA SUFICIENCIA DE RANGO PARA p=",p](*;Break[]*),Print["Valores proprios matriz Ke=",Valores];Print["NO tenemos la suficiencia de rango para p=",p]]
];