5 #include "Epetra_ConfigDefs.h" 8 #include "Epetra_MpiComm.h" 10 #include "Epetra_SerialComm.h" 13 #include "Teuchos_CommandLineProcessor.hpp" 14 #include "Teuchos_StandardCatchMacros.hpp" 15 #include "Teuchos_ParameterList.hpp" 16 #include "Teuchos_XMLParameterListCoreHelpers.hpp" 19 #include <boost/math/special_functions/gamma.hpp> 23 int main(
int argc,
char *argv[]){
25 std::string xmlInFileName =
"";
26 std::string extraXmlFile =
"";
27 std::string xmlOutFileName =
"paramList.out";
29 Teuchos::CommandLineProcessor clp(
false);
30 clp.setOption(
"xml-in-file",&xmlInFileName,
"The XML file to read into a parameter list");
31 clp.setDocString(
"TO DO.");
33 Teuchos::CommandLineProcessor::EParseCommandLineReturn
34 parse_return = clp.parse(argc,argv);
35 if( parse_return != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL ) {
36 std::cout <<
"\nEnd Result: TEST FAILED" << std::endl;
41 MPI_Init(&argc, &argv);
42 Epetra_MpiComm Comm(MPI_COMM_WORLD);
44 Epetra_SerialComm Comm;
47 Teuchos::RCP<Teuchos::ParameterList> paramList = Teuchos::rcp(
new Teuchos::ParameterList);
48 if(xmlInFileName.length()) {
49 Teuchos::updateParametersFromXmlFile(xmlInFileName, inoutArg(*paramList));
53 paramList->print(std::cout,2,
true,
true);
56 Teuchos::RCP<neumannInnerSurface_StochasticPolyconvexHGO> interface
59 std::ifstream parameters_file_1, parameters_file_2, parameters_file_3, parameters_file_4;
61 std::string path1 =
"/Users/Brian/Documents/Thesis/Trilinos_results/arteries/gmrf_neumann/a100_gamma3_delta010/";
62 parameters_file_1.open(path1+
"w1.txt");
63 parameters_file_2.open(path1+
"w2.txt");
64 parameters_file_3.open(path1+
"w3.txt");
65 parameters_file_4.open(path1+
"w4.txt");
67 unsigned int n_cells_p1_med = 297828;
68 unsigned int n_nodes_p1_med = 58464;
70 std::string path =
"/Users/Brian/Documents/Thesis/Trilinos/arteries/mesh/connectivity_p1_media.txt";
71 interface->get_media(n_cells_p1_med,n_nodes_p1_med,path);
73 if (parameters_file_1.is_open() && parameters_file_2.is_open() && parameters_file_3.is_open() && parameters_file_4.is_open()){
76 for (
int i=0;
i<n_nodes_p1_med; ++
i){
77 parameters_file_1 >> interface->w1_gmrf(
i);
78 parameters_file_2 >> interface->w2_gmrf(
i);
79 parameters_file_3 >> interface->w3_gmrf(
i);
80 parameters_file_4 >> interface->w4_gmrf(
i);
84 parameters_file_1.close();
85 parameters_file_2.close();
86 parameters_file_3.close();
87 parameters_file_4.close();
90 std::cout <<
"Couldn't open one of the parameters_file.\n";
94 int n_local_cells = interface->Mesh->n_local_cells;
95 int n_gauss_cells = interface->Mesh->n_gauss_cells;
96 std::vector<int> local_gauss_points(n_local_cells*n_gauss_cells);
97 for (
unsigned int e_lid=0; e_lid<n_local_cells; ++e_lid){
98 e_gid = interface->Mesh->local_cells[e_lid];
99 for (
unsigned int gp=0; gp<n_gauss_cells; ++gp){
100 local_gauss_points[e_lid*n_gauss_cells+gp] = e_gid*n_gauss_cells+gp;
104 Epetra_Map GaussMap(-1,n_local_cells*n_gauss_cells,&local_gauss_points[0],0,Comm);
106 Epetra_Vector mu1_gmrf(GaussMap);
107 Epetra_Vector mu2_gmrf(GaussMap);
108 Epetra_Vector mu3_gmrf(GaussMap);
109 Epetra_Vector mu4_gmrf(GaussMap);
110 Epetra_Vector x_coord(GaussMap);
111 Epetra_Vector y_coord(GaussMap);
112 Epetra_Vector z_coord(GaussMap);
115 double xi,
eta, zeta;
116 Epetra_SerialDenseMatrix matrix_X(3,interface->Mesh->el_type);
117 Epetra_SerialDenseVector vector_X(3);
118 Epetra_SerialDenseVector N(interface->Mesh->el_type);
119 for (
unsigned int e_lid=0; e_lid<n_local_cells; ++e_lid){
120 e_gid = interface->Mesh->local_cells[e_lid];
121 for (
int inode=0; inode<interface->Mesh->el_type; ++inode){
122 node = interface->Mesh->cells_nodes[interface->Mesh->el_type*e_gid+inode];
123 matrix_X(0,inode) = interface->Mesh->nodes_coord[3*node+0];
124 matrix_X(1,inode) = interface->Mesh->nodes_coord[3*node+1];
125 matrix_X(2,inode) = interface->Mesh->nodes_coord[3*node+2];
127 for (
unsigned int gp=0; gp<n_gauss_cells; ++gp){
128 xi = interface->Mesh->xi_cells[gp];
129 eta = interface->Mesh->eta_cells[gp];
130 zeta = interface->Mesh->zeta_cells[gp];
132 vector_X.Multiply(
'N',
'N',1.0,matrix_X,N,0.0);
133 interface->get_material_parameters(e_lid,gp);
134 mu1_gmrf[GaussMap.LID(
int(e_gid*n_gauss_cells+gp))] = interface->mu1;
135 mu2_gmrf[GaussMap.LID(
int(e_gid*n_gauss_cells+gp))] = interface->mu2;
136 mu3_gmrf[GaussMap.LID(
int(e_gid*n_gauss_cells+gp))] = interface->mu3;
137 mu4_gmrf[GaussMap.LID(
int(e_gid*n_gauss_cells+gp))] = interface->mu4;
138 x_coord[GaussMap.LID(
int(e_gid*n_gauss_cells+gp))] = vector_X(0);
139 y_coord[GaussMap.LID(
int(e_gid*n_gauss_cells+gp))] = vector_X(1);
140 z_coord[GaussMap.LID(
int(e_gid*n_gauss_cells+gp))] = vector_X(2);
145 int NumTargetElements = 0;
146 if (Comm.MyPID()==0){
147 NumTargetElements = interface->Mesh->n_cells*n_gauss_cells;
149 Epetra_Map MapOnRoot(-1,NumTargetElements,0,Comm);
150 Epetra_Export ExportOnRoot(GaussMap,MapOnRoot);
151 Epetra_MultiVector lhs_root(MapOnRoot,
true);
152 lhs_root.Export(mu1_gmrf,ExportOnRoot,Insert);
153 std::string
filename =
"/Users/Brian/Documents/Thesis/Trilinos_results/arteries/boost_gmrf/gmrf_gauss/mu1_gmrf.mtx";
154 error = EpetraExt::MultiVectorToMatrixMarketFile(filename.c_str(),lhs_root,0,0,
false);
156 lhs_root.PutScalar(0.0);
157 lhs_root.Export(mu2_gmrf,ExportOnRoot,Insert);
158 filename =
"/Users/Brian/Documents/Thesis/Trilinos_results/arteries/boost_gmrf/gmrf_gauss/mu2_gmrf.mtx";
159 error = EpetraExt::MultiVectorToMatrixMarketFile(filename.c_str(),lhs_root,0,0,
false);
161 lhs_root.PutScalar(0.0);
162 lhs_root.Export(mu3_gmrf,ExportOnRoot,Insert);
163 filename =
"/Users/Brian/Documents/Thesis/Trilinos_results/arteries/boost_gmrf/gmrf_gauss/mu3_gmrf.mtx";
164 error = EpetraExt::MultiVectorToMatrixMarketFile(filename.c_str(),lhs_root,0,0,
false);
166 lhs_root.PutScalar(0.0);
167 lhs_root.Export(mu4_gmrf,ExportOnRoot,Insert);
168 filename =
"/Users/Brian/Documents/Thesis/Trilinos_results/arteries/boost_gmrf/gmrf_gauss/mu4_gmrf.mtx";
169 error = EpetraExt::MultiVectorToMatrixMarketFile(filename.c_str(),lhs_root,0,0,
false);
171 lhs_root.PutScalar(0.0);
172 lhs_root.Export(x_coord,ExportOnRoot,Insert);
173 filename =
"/Users/Brian/Documents/Thesis/Trilinos_results/arteries/boost_gmrf/gmrf_gauss/x_coord.mtx";
174 error = EpetraExt::MultiVectorToMatrixMarketFile(filename.c_str(),lhs_root,0,0,
false);
176 lhs_root.PutScalar(0.0);
177 lhs_root.Export(y_coord,ExportOnRoot,Insert);
178 filename =
"/Users/Brian/Documents/Thesis/Trilinos_results/arteries/boost_gmrf/gmrf_gauss/y_coord.mtx";
179 error = EpetraExt::MultiVectorToMatrixMarketFile(filename.c_str(),lhs_root,0,0,
false);
181 lhs_root.PutScalar(0.0);
182 lhs_root.Export(z_coord,ExportOnRoot,Insert);
183 filename =
"/Users/Brian/Documents/Thesis/Trilinos_results/arteries/boost_gmrf/gmrf_gauss/z_coord.mtx";
184 error = EpetraExt::MultiVectorToMatrixMarketFile(filename.c_str(),lhs_root,0,0,
false);
int main(int argc, char *argv[])
void shape_functions(Epetra_SerialDenseVector &N, double &xi, double &eta, double &zeta)